Map matching method for distinguishing urban road context environment based on gnss multipath error

Abstract

The application provides a map matching method based on GNSS multipath error to distinguish urban road context environment, relates to global navigation satellite system positioning technology field, and includes the following steps: constructing a virtual city and satellite environment, simulating GNSS signal propagation by using a ray tracing algorithm, and generating a GNSS multipath error distribution heat map; based on error statistical characteristics, the road context environment is divided into multiple error level regions; a buffer region is defined with the GNSS trajectory point to be matched as the center, and then all candidate road segments in which the current GNSS trajectory point of the vehicle is located are screened out; through weighted calculation on the candidate road segments, the road segment with the highest path score is finally selected as the best matching road segment. The application quantitatively measures the spatial distribution of urban multipath error in advance, enhances the adaptive ability of the algorithm to the complex road context environment, effectively suppresses the matching deviation caused by the positioning jump in the high error area, and improves the accuracy and robustness of the map matching.

Description

Technical Field

[0001] This invention relates to the field of global navigation satellite system positioning technology, and in particular to a map matching method based on GNSS multipath error to distinguish urban road context. Background Technology

[0002] With the widespread adoption of Global Navigation Satellite Systems (GNSS) and location-based services, vehicle trajectory map matching technology plays an indispensable role in intelligent transportation, autonomous driving, and traffic flow monitoring. However, in complex urban environments (such as commercial districts with towering buildings or urban canyons), GNSS signals are easily blocked, reflected, or refracted by surrounding buildings, resulting in significant multipath effects. This interference between indirect and direct signals can lead to positioning jumps of tens of meters and severe deviations, becoming a core bottleneck restricting high-precision urban map matching.

[0003] Traditional map matching algorithms (such as nearest neighbor matching based on geometric distance and basic topology matching) perform poorly when processing such high-noise GNSS data. When vehicles are traveling in complex road networks (such as areas with parallel auxiliary roads or multi-level intersections) and encounter multipath errors, traditional algorithms that rely solely on single geometric distance and simple driving continuity are easily misled by positioning drift, incorrectly matching vehicles to adjacent but non-connected road segments, and lack error correction capabilities.

[0004] To mitigate the impact of multipath errors, existing technologies typically employ filtering algorithms (such as Kalman filtering and particle filtering) or multi-sensor fusion (such as GNSS combined with an inertial navigation system (INS)). However, filtering algorithms exhibit significant hysteresis and prediction bias when handling highly dynamic trajectories and complex multipath abrupt changes; while inertial navigation systems can maintain accuracy for short periods, they face the problem of error accumulation over time, and the large-scale deployment of high-precision hardware is prohibitively expensive. Furthermore, some machine learning-based error correction methods often rely on massive amounts of training data specific to particular scenarios, resulting in extremely limited generalization ability and environmental adaptability.

[0005] In summary, most existing map matching and error suppression techniques passively smooth out existing GNSS errors or rely on expensive additional hardware to compensate for them, failing to proactively perceive and quantify the local distortions caused by urban 3D space to GNSS signals at the source. Therefore, there is an urgent need in this field for a map matching method that can actively distinguish urban road context and adaptively adjust matching weights based on multi-path error distribution characteristics, in order to overcome the bottlenecks in matching accuracy and robustness in complex urban scenarios from the algorithmic level. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention analyzes the impact of GNSS multipath error on positioning accuracy by constructing virtual city environments and virtual satellite environments, and designs a map matching method based on GNSS multipath error to distinguish urban road context environments, thereby improving the accuracy and robustness of map matching.

[0007] To achieve the above objectives, this invention provides a map matching method for distinguishing urban road context based on GNSS multipath error, comprising the following steps: S1. Obtain GNSS multipath error distribution: Obtain a 3D model of the city to construct a virtual city environment; obtain ephemeris data to construct a virtual satellite environment; use ray tracing algorithm to simulate the propagation of GNSS signals in the virtual city environment and calculate the pseudorange deviation caused by multipath effects; combine GNSS measurement noise characteristics and positioning error model to create a GNSS multipath error distribution heatmap that characterizes the intensity and distribution of multipath effects in the city environment. S2. Delineate urban road context based on GNSS multipath error distribution heatmap: Extract statistical features such as mean, variance, skewness, kurtosis and distribution characteristics of GNSS multipath errors from the GNSS multipath error distribution heatmap; divide the urban road context into three categories of multipath error level regions: low, medium and high, based on the statistical features of GNSS multipath errors. S3. Determine candidate road segments based on geometric and error environment: Delineate a buffer area centered on the current GNSS trajectory point P to be matched; query the road network data within the buffer area. If a road segment is located within the buffer area and has spatial overlap with any multipath error level area, then determine it as the initial candidate road segment set for the current GNSS trajectory point. S4. Filter candidate road segments based on topology and driving logic: Further filter the initial candidate road segment set of the current GNSS trajectory points based on the topology of the road network and the driving logic of the vehicle to eliminate unreasonable candidate road segments. S5. Dynamically adjust weight parameters based on road context: dynamically adjust weight parameters according to the multi-path error level region where each candidate road segment is located; S6. Calculate the state probability: Based on GNSS trajectory data, road network data and dynamically adjusted weight parameters, calculate the observation probability of the current GNSS trajectory point on the candidate road segment and the transition probability between candidate road segments; S7. Solving the optimal path based on dynamic programming algorithm: The map matching problem is constructed as a hidden Markov model. The GNSS trajectory point sequence is used as the observation sequence and the candidate road segment sequence is used as the hidden state sequence. The state transition matrix is ​​established by combining the observation probability and the transition probability. The dynamic programming algorithm is used to find the optimal path and output the candidate road segment sequence corresponding to the path with the highest probability as the optimal matching road sequence for the vehicle.

[0008] Optionally, step S1 specifically includes: Step S101: Obtain the virtual 3D environment of the city By collecting urban 3D terrain and building models through oblique photography and LiDAR point cloud acquisition, a fine-grained virtual 3D urban environment is constructed. Step S102: Construct a virtual satellite environment Obtain ephemeris data from a GNSS repository or an actual GNSS receiver; The satellite's position in the geocentric fixed coordinate system is calculated using orbital parameters from ephemeris data; For each satellite, the ray tracing algorithm simulates two rays originating from the GNSS receiver: one is the line-of-sight (LOS) direct signal connecting the satellite and the receiver, and the other is the direction of the reflected signal. By calculating the length of these line segments, the existence of multipath error is detected and the additional distance of the reflected signal is calculated, thereby simulating the propagation path of the GNSS signal in a virtual city environment. Step S103: Calculation of multipath pseudorange deviation Based on known observation conditions and application scenarios, the additional path difference of multipath signals compared to line-of-sight (LOS) direct signals... The following three progressive models are used for calculation: Benchmark calculations based on physical definitions: , Local geometric calculations based on receiver perspective: , Engineering application calculations based on satellite elevation angles: , in, This represents the straight-line spatial reflection distance from the reflection point on the building's surface to the receiver. This represents the additional incident distance from the equiphase surface of the signal to the building's reflection point. This represents the geometric angle between the reflected path at the receiver and the equiphase surface of the signal, used to calculate the vertical component of the path. This represents the satellite's elevation angle, used to calculate the geometric relationship of the reflection path; The stochastic process of correlated GNSS noise is simulated using the Ornstein-Uhlenbeck process, and its stochastic differential equation is defined as follows: , in, and The parameters can be adjusted to match the randomness of a specific GNSS receiver model. The mean of the Ornstein-Uhlenbeck process. This represents the GNSS noise status value. For time differentiation, For the increment of the Wiener process (Brownian motion); Step S104: Generation of multipath error heatmap The positioning process involves calculating pseudorange measurements of the receiver and using these measurements to estimate the receiver's position. Evaluate the positioning accuracy of GNSS receivers in urban environments, including statistical modeling of errors, calculation of DOP values, and evaluation of positioning solution quality; Taking into account measurement noise, positioning, and positioning accuracy, multiple sampling points are selected within the target area. The GNSS positioning error intensity of each sampling point is calculated. The urban area is divided into grids, and the inverse distance weighted interpolation method is used to interpolate the GNSS positioning error intensity of the sampling points to the entire area to create a heat map to highlight and quantify the multipath effect in the urban environment.

[0009] Optionally, step S2 specifically includes: Step S201: Extract statistical features of GNSS multipath errors based on the GNSS multipath error distribution heatmap. Clean the relevant data of the GNSS multipath error distribution heatmap to remove obvious outliers and erroneous data; The area covered by the GNSS multipath error distribution heatmap is divided into grids of a certain size. The corresponding GNSS multipath error values ​​are extracted from each grid, and the mean, variance, skewness, and kurtosis of the GNSS multipath error are calculated. Step S202: Differentiate urban road context based on multipath error intensity K-means or Gaussian mixture models were used to cluster the mean, variance, skewness, and kurtosis of GNSS multipath errors. The statistical feature thresholds are dynamically adjusted based on the clustering results. The urban road context is then divided into three regions based on these thresholds: a high error intensity region, a medium error intensity region, and a low error intensity region.

[0010] Optionally, step S4 specifically includes: S401, Road Connectivity Filtering Based on the topology of the road network, candidate road segments connected to the road segment where the current GNSS trajectory point is located are selected. By analyzing the connectivity of the roads, candidate road segments not connected to the road segment where the current GNSS trajectory point is located are excluded. S402, Verification of the Reasonableness of Driving By combining vehicle speed, timestamp information, and historically matched locations, it is determined whether traveling from the previous location to the current candidate road segment conforms to the vehicle's driving logic.

[0011] Optionally, step S5 specifically includes: S501, Calculation of Relative Position Weights , , in, Indicates the instantaneous driving direction of the vehicle relative to the direction of traffic on the candidate road section The included angle, Dynamic relative position coefficients tied to multipath error levels, including relative position coefficients in the low error zone. Relative position coefficient of the error zone Relative position coefficient in high error zone And satisfy ; S502, Calculation of Weights for Correlation Between Previous and Previous Stages , in, This indicates the distance between the current candidate road segment and the historically matched road segments. Dynamic correlation coefficients linked to multipath error levels, including correlation coefficients in the low error region. Correlation coefficient before and after the mean error zone Correlation coefficient before and after the high error region And satisfy .

[0012] Optionally, the dynamic relative position coefficient Dynamic correlation coefficient Set it in any of the following ways: The error is obtained through regression analysis based on the matching error between historical trajectory data and the actual path; The parameters were optimized in different urban areas using grid search or Bayesian optimization methods. The settings are based on empirical values.

[0013] Optionally, step S6 specifically includes: Step S601: Calculate the observation probability Calculate the actual location of the vehicle in the candidate road segment GNSS trajectory points observed at the time The likelihood is calculated using the following formula: , in, The vertical distance distribution function is... GNSS trajectory points To candidate road sections vertical distance, Based on candidate road segments The observation confidence coefficient at the multipath error level, including the confidence coefficient of low-error observations. , Observation reliability coefficient High error observation reliability coefficient And satisfy ; Step S602: Calculate the transition probability Calculate the vehicle's journey from the previous time segment. Driving to the candidate road segment at the current moment The probability of is given by the following formula: , in, road_class This is the basic transition probability function calculated based on the connectivity and road grade between candidate road segments and historically matched road segments.

[0014] Optionally, the dynamic programming algorithm in step S7 is the Viterbi algorithm; Accordingly, a dynamic programming algorithm is used to optimize the path, and the sequence of candidate road segments corresponding to the path with the highest probability is output as the optimal matching road sequence for vehicle travel, including: Step S701, Initialization For the first trajectory point, initialize the initial path probability for each candidate road segment. and backtrack pointer Let the first trajectory point have Given 10 candidate road segments, calculate the initial path probability for each candidate road segment: , in, Initially in the candidate road segment And the first trajectory point was observed. The probability, Represents the first trajectory point's... One candidate road segment, ; Observation probability Initialized to 0, backtrack pointer Initialize to 0; Step S702, Recursive Calculation For each subsequent trajectory point, iterate through each of its candidate road segments. Calculate the arrival time of candidate road segments from the initial time. The maximum probability path, traversing the th All candidate road segments of each trajectory point (Right now The recursive formula is: , , in, For the first The trajectory points reach the candidate road segment The maximum probability among all possible paths, The indexes of candidate road segments that could generate this maximum transition probability are recorded. , , It is the first The total number of candidate road segments for each trajectory point; Step S703, Terminate Operation In all After processing all trajectory points, the road segment with the highest probability at the final moment is found and used as the endpoint of the optimal path:

[0015] in, For the last trajectory point The total number of candidate road segments; Step S704, Path Backtracking Initialize the optimal path sequence ; Based on the recorded backtracking pointer Backtracking from the endpoint: , The path scoring function of the Viterbi algorithm is used to calculate the score of the candidate road segment at each time step: , in, The length of the trajectory point sequence. This represents the t-th trajectory point. This represents the candidate road segment corresponding to the t-th trajectory point; Determine the optimal candidate road segment corresponding to each trajectory point at each time step, and finally output the sequence. This is the globally optimal vehicle map matching road sequence.

[0016] By adopting the above technical solution, the present invention has at least the following beneficial effects: This invention innovatively utilizes a virtual 3D urban model to accurately acquire the GNSS positioning error intensity within a target area in a physical simulation environment. The obtained data is then presented visually as a heatmap, clearly demonstrating the distribution of GNSS positioning errors. Based on this error distribution, a corresponding buffer zone is established, and all candidate road segments within this buffer zone are selected based on the vehicle's current GNSS trajectory point. Subsequently, by weighting these candidate road segments, the road segment with the highest path score is selected as the optimal matching road segment. This process enables the calculation of multipath effects in the urban environment before map matching, analyzing the impact area of ​​multipath errors, distinguishing the urban road context, and thus improving the map matching algorithm. The optimal matching road is obtained using weighted calculations to mitigate the impact of multipath errors on map matching. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 A flowchart of a map matching method for distinguishing urban road context based on GNSS multipath error is provided in an embodiment of the present invention. Figure 2 A schematic diagram illustrating the principle of GNSS signal ray tracing in an urban environment; Figure 3 A schematic diagram of the benchmark computational geometry based on physical definitions; Figure 4 This is a schematic diagram of the geometric principle of local geometric computation based on the receiver's perspective; Figure 5 A schematic diagram illustrating the computational geometry principle for engineering applications based on satellite elevation angles; Figure 6 This is a schematic diagram of the initial spatial overlay of candidate road segments based on the GNSS multipath error distribution heatmap provided in an embodiment of the present invention; Figure 7 This is a schematic diagram of the candidate road segment topology and probability transition of a Hidden Markov Model (HMM) provided in an embodiment of the present invention. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] like Figure 1 As shown, this invention provides a map matching method for distinguishing urban road context based on GNSS multipath error, comprising the following steps: S1. Obtain GNSS multipath error distribution In this embodiment, a virtual urban environment is constructed by acquiring a 3D model of the city; a virtual satellite environment is constructed by acquiring ephemeris data; the propagation of GNSS signals in the virtual urban environment is simulated using a ray tracing algorithm, and the pseudorange deviation caused by multipath effects is calculated; and a GNSS multipath error distribution heatmap that characterizes the intensity and distribution of multipath effects in the urban environment is created by combining GNSS measurement noise characteristics and positioning error models.

[0021] Specifically, the following steps are included: Step S101: Obtain the virtual 3D environment of the city By collecting urban 3D terrain and building models through oblique photography and LiDAR point cloud acquisition, a fine-grained urban virtual 3D environment is constructed, with the model resolution preferably not less than 0.5 meters.

[0022] Step S102: Construct a virtual satellite environment

[0023] Ephemeris data is obtained from a GNSS repository or an actual GNSS receiver. Ephemeris data contains Kepler features of the satellite orbit, such as orbital inclination and eccentricity.

[0024] The satellite's position in the geocentric fixed coordinate system is calculated using orbital parameters from ephemeris data.

[0025] For each satellite, the ray tracing algorithm simulates two rays originating from the GNSS receiver: one is the line-of-sight (LOS) direct signal connecting the satellite and the receiver, and the other is the direction of the reflected signal. By calculating the lengths of these line segments, the existence of multipath errors is detected, and the additional distance of the reflected signal is calculated, thereby simulating the propagation path of the GNSS signal in a virtual city environment.

[0026] Step S103: Calculation of multipath pseudorange deviation

[0027] The main physical effect of multipath signals is that it increases the number of propagation paths of the signal, such as... Figure 2As shown. Based on known observation conditions and application scenarios, the additional path difference of multipath signals compared to line-of-sight (LOS) direct signals. The following three progressive models are used for calculation: Assuming the GNSS satellite is extremely far from Earth, the line-of-sight (LOS) direct signal reaching the observation area and the incident signal causing multipath propagation can be considered as parallel beams. A perpendicular line is drawn from the receiver towards the direction of the incident signal; this perpendicular line physically represents the equiphase surface of the signal.

[0028] Benchmark calculations based on physical definitions: such as Figure 3 As shown, when the additional incident distance from the equiphase surface of the signal to the building reflection point is theoretically known... When, the formula is as follows: , Local geometric calculation based on receiver perspective: when it cannot be directly measured However, the geometric angle between the reflection path at the receiver and the equiphase surface of the signal is known. At times, such as Figure 4 As shown, in the right triangle formed by the equiphase surface, the incident path, and the reflected path, according to the sine theorem, we know... Substituting into the formula, we get: , Engineering applications based on satellite elevation angles: In practical engineering calculations, satellite elevation angles can only be obtained through ephemeris data. At times, such as Figure 5 This means that the included angle can be derived using the angular relationships of plane geometry. With satellite elevation angle The conversion relationship is According to the Law of Sines, we know Substituting this relationship into the formula, we get: , in, This represents the straight-line spatial reflection distance from the reflection point on the building's surface to the receiver. This represents the additional incident distance from the equiphase surface of the signal to the building's reflection point. This represents the geometric angle between the reflected path at the receiver and the equiphase surface of the signal, used to calculate the vertical component of the path. This represents the satellite's elevation angle, used to calculate the geometric relationship of the reflection path.

[0029] The stochastic process of correlated GNSS noise is simulated using the Ornstein-Uhlenbeck process, and its stochastic differential equation is defined as follows: , in, and The parameters can be adjusted to match the randomness of a specific GNSS receiver model. The mean of the Ornstein-Uhlenbeck process. This represents the GNSS noise status value. For time differentiation, This is the increment of the Wiener process (Brownian motion).

[0030] Step S104: Generation of multipath error heatmap

[0031] The positioning process involves calculating pseudorange measurements of the receiver and using these measurements to estimate the receiver's position.

[0032] Evaluating the positioning accuracy of GNSS receivers in urban environments includes statistical modeling of errors, calculation of DOP values, and assessment of positioning solution quality.

[0033] Taking into account measurement noise, positioning, and positioning accuracy, multiple sampling points are selected within the target area. The GNSS positioning error intensity of each sampling point is calculated. The urban area is divided into grids, and the inverse distance weighted interpolation method is used to interpolate the GNSS positioning error intensity of the sampling points to the entire area to create a heat map to highlight and quantify the multipath effect in the urban environment.

[0034] S2. Defining the urban road context environment based on GNSS multipath error distribution heatmap

[0035] In this embodiment, statistical features such as mean, variance, skewness, kurtosis and distribution characteristics of GNSS multipath errors are extracted from the GNSS multipath error distribution heatmap; based on the statistical features of GNSS multipath errors, the urban road context environment is divided into three categories of multipath error level regions: low, medium and high.

[0036] Step S2 specifically includes: Step S201: Extract statistical features of GNSS multipath errors based on the GNSS multipath error distribution heatmap. Clean the GNSS multipath error distribution heatmap data to remove obvious outliers and erroneous data.

[0037] The area covered by the GNSS multipath error distribution heatmap is divided into grids of a certain size. The corresponding GNSS multipath error value is extracted from each grid. Different colors in the heatmap represent different error magnitudes. These values ​​are then converted into specific numerical forms to calculate the mean, variance, skewness, and kurtosis of the GNSS multipath error. The mean reflects the central tendency of the multipath error in the region, the variance measures the dispersion of the multipath error, i.e., the degree of deviation between the multipath error value and the mean, the skewness describes the symmetry of the multipath error distribution, and the kurtosis reflects the peak of the error distribution.

[0038] Step S202: Differentiate urban road context based on multipath error intensity

[0039] K-means or Gaussian mixture models were used to cluster the mean, variance, skewness, and kurtosis of GNSS multipath errors.

[0040] Based on the clustering results, the statistical feature thresholds are dynamically adjusted, and the urban road context is divided into three intervals according to the statistical feature thresholds: mean error: <5m is the low error zone, 5~15m is the medium error zone, and >15m is the high error zone; variance: <2m² is the low error zone, 2~8m² is the medium error zone, and >8m² is the high error zone; skewness and kurtosis are used to help judge the concentration and extreme nature of the error distribution. Finally, the area is divided into high error intensity area (red area), medium error intensity area (yellow area), and low error intensity area (green area).

[0041] S3. Determining candidate road segments based on geometric and error environment

[0042] In this embodiment, as Figure 6 As shown, a buffer zone is defined with the current GNSS trajectory point P to be matched as the center; the road network data within the buffer zone is queried. If a road segment is located within the buffer zone and has spatial overlap with any multipath error level zone, it is determined as the initial candidate road segment set for the current GNSS trajectory point.

[0043] The buffer zone can be defined using any of the following methods: Fixed radius buffer: A buffer with a fixed radius is constructed, centered on the GNSS trajectory point and with an error threshold as the radius. All road segments within the buffer are considered candidate road segments and participate in subsequent map matching calculations.

[0044] Dynamic buffer adjustment: Considering that the error distribution may differ in different areas, a dynamic buffer adjustment method can be adopted. The initial buffer radius is set to the error threshold. If the initial candidate road segment set is empty, the radius is increased by a fixed step size until the initial candidate road segment set contains at least one road segment.

[0045] S4. Filtering candidate road segments based on topology and driving logic

[0046] In this embodiment, the initial candidate road segment set of the current GNSS trajectory point is further filtered based on the topology of the road network and the driving logic of the vehicle to eliminate unreasonable candidate road segments.

[0047] Step S4 specifically includes: S401, Road Connectivity Filtering Based on the topology of the road network, candidate road segments connected to the road segment where the current GNSS trajectory point is located are selected. By analyzing the connectivity of the roads, candidate road segments that are not connected to the road segment where the current GNSS trajectory point is located are eliminated, thereby reducing the number of candidate road segments.

[0048] S402, Verification of the Reasonableness of Driving

[0049] By combining vehicle speed, timestamp information, and historically matched locations, it is determined whether traveling from the previous location to the current candidate road segment conforms to the vehicle's driving logic.

[0050] S5. Dynamically adjust weight parameters based on road context.

[0051] In this embodiment, the weight parameters are dynamically adjusted according to the multipath error level region where each candidate road segment is located.

[0052] Step S5 specifically includes: S501. The relative position weight is calculated based on the angle between the vehicle's instantaneous driving direction and the candidate road segment's travel direction, and a dynamic relative position coefficient bound to the error level is introduced to achieve the calculation. The formula is as follows: , , in, Indicates the instantaneous driving direction of the vehicle relative to the direction of traffic on the candidate road section The included angle, Dynamic relative position coefficients tied to multipath error levels, including relative position coefficients in the low error zone. Relative position coefficient of the error zone Relative position coefficient in high error zone And satisfy .

[0053] The smaller the angle between the vehicle's instantaneous driving direction and the travel direction of each candidate road segment, the greater the weight of the factor representing the relative position between the positioning point and the candidate road segment, and the greater the probability that the corresponding road segment is the best matching road.

[0054] S502. Calculation is based on the spatial distance between the current candidate road segment and historically matched successful road segments, and a dynamic correlation coefficient bound to the error level is introduced to realize the correlation weight calculation, as shown in the following formula: , in, This indicates the distance between the current candidate road segment and the historically matched road segments. Dynamic correlation coefficients linked to multipath error levels, including correlation coefficients in the low error region. Correlation coefficient before and after the mean error zone Correlation coefficient before and after the high error region And satisfy .

[0055] The correlation weight value reflects the degree of correlation between the current candidate road segment and the historical matching road segments. The higher the degree of correlation, the greater the weight value.

[0056] Dynamic relative position coefficient Dynamic correlation coefficient The parameters can be set using any of the following methods: regression analysis based on the matching error between historical trajectory data and the actual path; parameter tuning in different urban areas using grid search or Bayesian optimization methods; or empirical values. =0.5, =1.0, =2.0; =1.5, =1.0, =0.5.

[0057] S6. Calculate the state probability

[0058] In this embodiment, based on GNSS trajectory data, road network data, and dynamically adjusted weight parameters, combined with... Figure 4 The diagram illustrates the topology and probability transition relationships of candidate road segments in a Hidden Markov Model (HMM). Two core probabilities in the HMM are calculated: the observation probability of the current GNSS trajectory point on a candidate road segment and the transition probability between candidate road segments. These probabilities will be used in subsequent dynamic programming algorithms to solve for the optimal matching path.

[0059] Step S6 specifically includes: Step S601: Based on the vertical distance distribution function from the current GNSS trajectory point to the candidate road segment, and introducing the observation confidence coefficient based on the error level and the relative position weight as adjustment factors, calculate the distance when the vehicle is actually located on the candidate road segment. GNSS trajectory points observed at the time The likelihood is calculated using the following formula: , in, The vertical distance distribution function is... GNSS trajectory points To candidate road sections vertical distance, Based on candidate road segments The observation confidence coefficient at the multipath error level, including the confidence coefficient of low-error observations. , Observation reliability coefficient High error observation reliability coefficient And satisfy .

[0060] The calculation of this probability takes into account both geometric relationships and the current road context. The calculation of the observation probability consists of three parts: the basic observation likelihood. Describes the degree of geometric matching between trajectory points and road segments under ideal conditions. Relative position weights. The importance of dynamically adjusting orientation matching. Observational reliability coefficient. Based on the multipath error level, the credibility of the current GNSS trajectory point itself is globally modulated.

[0061] Basic observation likelihood function A typical implementation is as follows: , in, This is an adjustable parameter used to control the sensitivity to distance. Since HMM solutions typically rely on relative probabilities, in practical calculations, it is necessary to apply this parameter to the same trajectory points. The observed likelihoods of all candidate road segments are normalized to form a probability distribution: , in, For trajectory points The total number of candidate road segments.

[0062] Step S602: Calculate the transition probability

[0063] Based on the fundamental transition probability function between candidate road segments and historically matched road segments, and introducing preceding and following correlation weights as adjustment factors, the vehicle's transition from the road segment at the previous time step is calculated. Driving to the candidate road segment at the current moment The probability of is given by the following formula: , in, road_class This is the basic transition probability function calculated based on the connectivity and road grade between candidate road segments and historically matched road segments.

[0064] This probability reflects the continuity and rationality of the path in time and space.

[0065] Basic transition probability function The rules are defined as follows: Connectivity determination: If the road segment and If not directly connected, then =0.

[0066] Road grade bonus: If two road sections are connected, then... The value of is positively correlated with the road grade.

[0067] , in, >0.

[0068] Similar to observation probability, for a given predecessor state It is necessary to consider all possible successor states. The transition likelihood is normalized to form a probability distribution: .

[0069] S7. Solving the optimal path based on dynamic programming algorithm

[0070] In this embodiment, as Figure 7 As shown, the map matching problem is constructed as a Hidden Markov Model, with the GNSS trajectory point sequence as the observation sequence and the candidate road segment sequence as the hidden state sequence. The state transition matrix is ​​established by combining the observation probability and the transition probability. The Viterbi algorithm is used to find the hidden state sequence most likely to generate the current observation sequence, that is, the optimal matching road sequence for vehicle travel.

[0071] Step S7 specifically includes: Step S701, Initialization For the first trajectory point, initialize the initial path probability for each candidate road segment. and backtrack pointer Let the first trajectory point have Given 10 candidate road segments, calculate the initial path probability for each candidate road segment: , in, Initially in the candidate road segment And the first trajectory point was observed. The probability, Represents the first trajectory point's... One candidate road segment, ; Since there are no prior initial state preferences and no predecessor states, the observation probability is directly expressed as... Initialize to 0.

[0072] Since this is the initial trajectory point, there is no predecessor state, so the backtracking pointer is directly set... Initialize to 0.

[0073] Step S702, Recursive Calculation

[0074] For each subsequent trajectory point, iterate through each of its candidate road segments. Calculate the arrival time of candidate road segments from the initial time. The maximum probability path, traversing the th All candidate road segments of each trajectory point (Right now The recursive formula is: , , in, For the first The trajectory points reach the candidate road segment The maximum probability among all possible paths, The indexes of candidate road segments that could generate this maximum transition probability are recorded. , , It is the first The total number of candidate road segments for each trajectory point; Step S703, Terminate Operation In all After processing all trajectory points, the road segment with the highest probability at the final moment is found and used as the endpoint of the optimal path:

[0075] in, For the last trajectory point The total number of candidate road segments; Step S704, Path Backtracking Initialize the optimal path sequence ; Based on the recorded backtracking pointer Backtracking from the endpoint: , The path scoring function of the Viterbi algorithm is used to calculate the score of the candidate road segment at each time step: , in, The length of the trajectory point sequence. This represents the t-th trajectory point. This represents the candidate road segment corresponding to the t-th trajectory point; Determine the optimal candidate road segment corresponding to each trajectory point at each time step, and finally output the sequence. This is the globally optimal vehicle map matching road sequence.

[0076] The present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

Claims

1. A map matching method for distinguishing urban road context based on GNSS multipath error, characterized in that, The steps include the following: S1. Obtain GNSS multipath error distribution: Obtain a 3D model of the city to construct a virtual city environment; obtain ephemeris data to construct a virtual satellite environment; use ray tracing algorithm to simulate the propagation of GNSS signals in the virtual city environment and calculate the pseudorange deviation caused by multipath effects; combine GNSS measurement noise characteristics and positioning error model to create a GNSS multipath error distribution heatmap that characterizes the intensity and distribution of multipath effects in the city environment. S2. Delineate the urban road context based on the GNSS multipath error distribution heatmap: Extract statistical features such as mean, variance, skewness, kurtosis and distribution characteristics of GNSS multipath errors from the GNSS multipath error distribution heatmap; Based on the statistical characteristics of GNSS multipath error, the urban road context environment is divided into three categories of multipath error levels: low, medium, and high. S3. Determine candidate road segments based on geometric and error environment: Delineate a buffer area centered on the current GNSS trajectory point P to be matched; query the road network data within the buffer area. If a road segment is located within the buffer area and has spatial overlap with any multipath error level area, then determine it as the initial candidate road segment set for the current GNSS trajectory point. S4. Filter candidate road segments based on topology and driving logic: Further filter the initial candidate road segment set of the current GNSS trajectory points based on the topology of the road network and the driving logic of the vehicle to eliminate unreasonable candidate road segments. S5. Dynamically adjust weight parameters based on road context: dynamically adjust weight parameters according to the multi-path error level region where each candidate road segment is located; S6. Calculate the state probability: Based on GNSS trajectory data, road network data and dynamically adjusted weight parameters, calculate the observation probability of the current GNSS trajectory point on the candidate road segment and the transition probability between candidate road segments; S7. Solving the optimal path based on dynamic programming algorithm: The map matching problem is constructed as a hidden Markov model. The GNSS trajectory point sequence is used as the observation sequence and the candidate road segment sequence is used as the hidden state sequence. The state transition matrix is ​​established by combining the observation probability and the transition probability. The dynamic programming algorithm is used to find the optimal path and output the candidate road segment sequence corresponding to the path with the highest probability as the optimal matching road sequence for the vehicle.

2. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 1, characterized in that, Step S1 specifically includes: Step S101: Obtain the virtual 3D environment of the city By collecting urban 3D terrain and building models through oblique photography and LiDAR point cloud acquisition, a fine-grained virtual 3D urban environment is constructed. Step S102: Construct a virtual satellite environment Obtain ephemeris data from a GNSS repository or an actual GNSS receiver; The satellite's position in the geocentric fixed coordinate system is calculated using orbital parameters from ephemeris data; For each satellite, the ray tracing algorithm simulates two rays originating from the GNSS receiver: one is the line-of-sight (LOS) direct signal connecting the satellite and the receiver, and the other is the direction of the reflected signal. By calculating the length of these line segments, the existence of multipath error is detected and the additional distance of the reflected signal is calculated, thereby simulating the propagation path of the GNSS signal in a virtual city environment. Step S103: Calculation of multipath pseudorange deviation Based on known observation conditions and application scenarios, the additional path difference of multipath signals compared to line-of-sight (LOS) direct signals... The following three progressive models are used for calculation: Benchmark calculations based on physical definitions: ， Local geometric calculations based on receiver perspective: ， Engineering application calculations based on satellite elevation angles: ， in, This represents the straight-line spatial reflection distance from the reflection point on the building's surface to the receiver. This represents the additional incident distance from the equiphase surface of the signal to the building's reflection point. This represents the geometric angle between the reflected path at the receiver and the equiphase surface of the signal, used to calculate the vertical component of the path. This represents the satellite's elevation angle, used to calculate the geometric relationship of the reflection path; The stochastic process of correlated GNSS noise is simulated using the Ornstein-Uhlenbeck process, and its stochastic differential equation is defined as follows: ， in, and The parameters can be adjusted to match the randomness of a specific GNSS receiver model. The mean of the Ornstein-Uhlenbeck process. This represents the GNSS noise status value. For time derivative, For the increment of the Wiener process (Brownian motion); Step S104: Generation of multipath error heatmap The positioning process involves calculating pseudorange measurements of the receiver and using these measurements to estimate the receiver's position. Evaluate the positioning accuracy of GNSS receivers in urban environments, including statistical modeling of errors, calculation of DOP values, and evaluation of positioning solution quality; Taking into account measurement noise, positioning, and positioning accuracy, multiple sampling points are selected within the target area. The GNSS positioning error intensity of each sampling point is calculated. The urban area is divided into grids, and the inverse distance weighted interpolation method is used to interpolate the GNSS positioning error intensity of the sampling points to the entire area to create a heat map to highlight and quantify the multipath effect in the urban environment.

3. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 1, characterized in that, Step S2 specifically includes: Step S201: Extract statistical features of GNSS multipath errors based on the GNSS multipath error distribution heatmap. Clean the relevant data of the GNSS multipath error distribution heatmap to remove obvious outliers and erroneous data; The area covered by the GNSS multipath error distribution heatmap is divided into grids of a certain size. The corresponding GNSS multipath error values ​​are extracted from each grid, and the mean, variance, skewness, and kurtosis of the GNSS multipath error are calculated. Step S202: Differentiate urban road context based on multipath error intensity K-means or Gaussian mixture models were used to cluster the mean, variance, skewness, and kurtosis of GNSS multipath errors. The statistical feature thresholds are dynamically adjusted based on the clustering results. The urban road context is then divided into three regions based on these thresholds: a high error intensity region, a medium error intensity region, and a low error intensity region.

4. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 1, characterized in that, Step S4 specifically includes: S401, Road Connectivity Filtering Based on the topology of the road network, candidate road segments connected to the road segment where the current GNSS trajectory point is located are selected. By analyzing the connectivity of the roads, candidate road segments not connected to the road segment where the current GNSS trajectory point is located are excluded. S402, Verification of the Reasonableness of Driving By combining vehicle speed, timestamp information, and historically matched locations, it is determined whether traveling from the previous location to the current candidate road segment conforms to the vehicle's driving logic.

5. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 1, characterized in that, Step S5 specifically includes: S501, Calculation of Relative Position Weights ， ， in, Indicates the instantaneous driving direction of the vehicle relative to the direction of traffic on the candidate road section The included angle, Dynamic relative position coefficients tied to multipath error levels, including relative position coefficients in the low error zone. Relative position coefficient of the error zone Relative position coefficient in high error zone And satisfy ; S502, Calculation of Weights for Corresponding Relationships Before and After ， in, This indicates the distance between the current candidate road segment and the historically matched road segments. Dynamic correlation coefficients linked to multipath error levels, including correlation coefficients in the low error region. Correlation coefficient before and after the mean error zone Correlation coefficient before and after the high error region And satisfy .

6. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 5, characterized in that, The dynamic relative position coefficient Dynamic correlation coefficient Set it in any of the following ways: The error is obtained through regression analysis based on the matching error between historical trajectory data and the actual path; The parameters were optimized in different urban areas using grid search or Bayesian optimization methods. The settings are based on empirical values.

7. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 5, characterized in that, Step S6 specifically includes: Step S601: Calculate the observation probability Calculate the actual location of the vehicle in the candidate road segment GNSS trajectory points observed at the time The likelihood is calculated using the following formula: ， in, The vertical distance distribution function is... GNSS trajectory points To candidate road sections vertical distance, Based on candidate road segments The observation confidence coefficient at the multipath error level, including the confidence coefficient of low-error observations. , Observation reliability coefficient High error observation reliability coefficient And satisfy ; Step S602: Calculate the transition probability Calculate the vehicle's journey from the previous time segment. Driving to the candidate road segment at the current moment The probability of is given by the following formula: ， in, road_class This is the basic transition probability function calculated based on the connectivity and road grade between candidate road segments and historically matched road segments.

8. The map matching method for distinguishing urban road context based on GNSS multipath error according to claim 7, characterized in that, The dynamic programming algorithm used in step S7 is the Viterbi algorithm. Accordingly, a dynamic programming algorithm is used to optimize the path, and the sequence of candidate road segments corresponding to the path with the highest probability is output as the optimal matching road sequence for vehicle travel, including: Step S701, Initialization For the first trajectory point, initialize the initial path probability for each candidate road segment. and backtrack pointer Let the first trajectory point have Given 10 candidate road segments, calculate the initial path probability for each candidate road segment: ， in, Initially in the candidate road segment And the first trajectory point was observed. The probability, Represents the first trajectory point's... One candidate road segment, ; Observation probability Initialized to 0, backtrack pointer Initialize to 0; Step S702, Recursive Calculation For each subsequent trajectory point, iterate through each of its candidate road segments. Calculate the arrival time of candidate road segments from the initial time. The maximum probability path, traversing the th All candidate road segments of each trajectory point (Right now The recursive formula is: ， ， in, For the first The trajectory points reach the candidate road segment The maximum probability among all possible paths, The indexes of candidate road segments that can generate this maximum transition probability are recorded. , , It is the first The total number of candidate road segments for each trajectory point; Step S703, Terminate Operation In all After processing all trajectory points, the road segment with the highest probability at the final moment is found and used as the endpoint of the optimal path: in, For the last trajectory point The total number of candidate road segments; Step S704, Path Backtracking Initialize the optimal path sequence ; Based on the recorded backtracking pointer Backtracking from the endpoint: ， The path scoring function of the Viterbi algorithm is used to calculate the score of the candidate road segment at each time step: ， in, The length of the trajectory point sequence. This represents the t-th trajectory point. This represents the candidate road segment corresponding to the t-th trajectory point; Determine the optimal candidate road segment corresponding to each trajectory point at each time step, and finally output the sequence. This is the globally optimal vehicle map matching road sequence.

Images