Low-cost positioning method and system for urban canyon environment based on 3D map

Abstract

The application provides a low-cost positioning method and system based on a 3D map of an urban canyon environment. The application uses high-precision two-dimensional image data provided by Google Earth to extract height information of high-rise buildings in the urban canyon area through a shadow height measurement method, and fuses the height information with building contour information obtained through OSM to construct a 3D map suitable for shadow matching positioning. Meanwhile, the shadow matching positioning algorithm is improved, the building information provided by the 3D map is used for shadow matching positioning calculation, and GNSS information is fused. The application can effectively overcome the limitation of the low popularization degree of high-precision 3D maps on the application of shadow matching positioning, improve the poor positioning accuracy of traditional shadow matching positioning in the cross-street direction, has no additional hardware cost, is simple to operate, can effectively avoid the storage and processing pressure of high-precision 3D maps on terminal devices, is suitable for low-cost terminal devices, and can effectively improve the positioning accuracy in the urban canyon environment.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of positioning, and relates to GNSS positioning technology, in particular to a low-cost positioning method and system for urban canyon environment based on a 3D map. BACKGROUND

[0002] GNSS(Global Navigation Satellite System,global navigation satellite system) positioning has become an indispensable part of people's daily life and industrial applications, the most common application is in vehicle navigation system, smart phone positioning service, map software, etc. In open areas, GNSS pseudorange single point positioning can usually reach meter-level accuracy, which can meet general positioning needs. But in urban canyon environment with high-rise buildings on both sides, the number of visible stars decreases sharply, and GNSS signals are easily blocked, reflected, etc. by high-rise buildings, resulting in signal attenuation, multipath interference, non-line-of-sight signal(Non Light Of Sight, NLOS) reception, and the resulting pseudorange single point positioning error can reach dozens of meters or even cannot be positioned. Introducing additional sensors such as laser radar, camera, IMU, etc. can improve positioning accuracy to some extent through multi-sensor fusion technology, but there are limitations in practical applications limited by cost, size, weight, etc. With the development of three-dimensional(three-dimensional, 3D) city maps, considering the geometric relationship between satellites and urban buildings, it provides a basis for judging LOS / NLOS signals to weaken multipath error in urban environment, and provides a new idea for improving GNSS positioning accuracy in urban canyon environment. The traditional shadow matching positioning algorithm uses 3D city map information, satellite ephemeris data and received satellite signal data for matching to narrow down the area where the receiving device may be located, thereby realizing the positioning of the receiving device, which can effectively improve the positioning accuracy in the direction of crossing the street in urban canyon environment, and has the advantages of low cost, small amount of calculation, etc., but the accuracy in the direction along the street is poor. In addition, the popularization degree of city 3D map with building height information is not high, and the method of establishing 3D map by using laser radar scanning, total station measurement, etc. has high cost and is complex, which is not conducive to the application of positioning method assisted by shadow matching algorithm in low-cost receiving device. Shadow matching positioning only needs to use the boundary information of the building, but the existing high-precision 3D city model often contains a large amount of additional detailed information, which increases the complexity of processing map data, and at the same time, storing a large amount of high-precision model data on low-cost terminal devices will bring additional pressure, affecting the running efficiency and resource utilization of the device. SUMMARY

[0003] To address the challenges of existing technologies, this paper provides a low-cost positioning method and system for urban canyon environments based on 3D maps. This method uses OpenStreetMap (OSM, an open source map project) and Google Earth to create a 3D map of the target area. Building information provided by the 3D map is then used to perform shadow matching positioning calculations to assist GNSS positioning. Furthermore, an improved method for the traditional shadow matching positioning algorithm is proposed to enhance the positioning accuracy of the shadow matching algorithm itself and integrate it with GNSS information to improve overall positioning accuracy.

[0004] In order to achieve the above object, the technical solution adopted by the present invention is:

[0005] A low-cost positioning method for urban canyon environments based on 3D maps includes the following steps:

[0006] Step 1: Extract building outline data:

[0007] Use map software to obtain the two-dimensional geographic data of points, lines and surfaces of the target area for which a 3D map needs to be created;

[0008] Step 2: Extract building height information data:

[0009] Obtain a two-dimensional image of the target area in step 1 through Google Earth, perform image processing on the two-dimensional image, measure the length of the building shadow in the processed two-dimensional image, and use shadow height measurement to estimate building height information based on the relationship between building height and building shadow length;

[0010] Step 3: Fusion of 2D building outlines and building height information:

[0011] Importing the building height information extracted in step 2 into the two-dimensional geographic data extracted in step 1 to generate a 3D map suitable for shadow matching positioning that includes the building height information;

[0012] Step 4: Determine the search area:

[0013] Use the terminal device to receive GNSS raw signal data, use pseudorange information and broadcast ephemeris to perform pseudorange single-point positioning, obtain the coordinates of the initial position P in the Earth-centered Earth-fixed coordinate system and convert them into Gaussian projection coordinates, determine the size of the search area radius r based on the quality of the satellite data, and generate a circular search area with P as the center and r as the radius;

[0014] Step 5: Generate candidate location grid:

[0015] Within the circular search area determined in step 4, candidate positions are set according to preset intervals to establish a candidate position grid;

[0016] Step 6: Generate a sky occlusion map for each candidate location:

[0017] Generate a sky occlusion map at each candidate location using building height information, the sky occlusion map containing the elevation angle information of the building boundary corresponding to each angle within a 360° azimuth range;

[0018] Step 7, predict satellite visibility:

[0019] Get satellite S through broadcast ephemeris i Azimuth and altitude angle Compare the jth candidate position C j Azimuth in the sky occlusion map Building height angle E B and satellite elevation angle The relationship between the predicted satellite S i At candidate position C j visibility of the place;

[0020] Step 8: Observe satellite visibility

[0021] Set a signal-to-noise ratio (SNR) or carrier-to-noise ratio (CNR) threshold based on empirical values. Determine the type of satellite signal and the visibility of the observed satellite by comparing the SNR or CNR of the received GNSS signal with the set threshold.

[0022] Step 9, matching scoring:

[0023] At each candidate location, score the match between the satellite visibility predicted in step 7 and the satellite visibility observed in step 8;

[0024] Step 10: Select the candidate position with the highest score to calculate the shadow matching positioning solution:

[0025] Count the total scores of all satellites at each candidate position after step 9, select the candidate position with the highest score, and take the average of the easting and northing coordinates of these candidate positions as the coordinates of the shadow matching positioning solution;

[0026] Step 11: Fuse the shadow matching positioning solution with the GNSS velocity information:

[0027] Based on the shadow matching positioning solution obtained in step 10, the GNSS velocity information is introduced, and the state equation and observation equation of the system are established through EKF design. The shadow matching positioning solution and GNSS velocity information are fused to obtain the optimal position estimate.

[0028] Furthermore, the point-line-surface two-dimensional geographic data of step 1 is obtained by the following process: selecting the target area for which a 3D map needs to be established in OpenStreetMap, importing the target area file into the map editing software ArcMap to extract the point-line-surface two-dimensional geographic data, wherein the point-line-surface two-dimensional geographic data includes road geographic information and building outlines.

[0029] Furthermore, step 2 includes the following sub-steps:

[0030] Step 2-1, running Google Earth, selecting the shooting time of an image in the historical images of Google Earth, selecting the target area range corresponding to step 1, and downloading the 2D image of the selected time and target area;

[0031] Step 2-2, performing image processing on the downloaded image;

[0032] Step 2-3, use the NOAA solar altitude angle online calculation tool to complete the calculation of the solar altitude angle β;

[0033] Steps 2-4: Select a typical high-rise building with a known height of H as a reference. Use the distance and area measurement tools in Google Earth to measure the length of the high-rise building's shadow as L. Based on the principle of high-rise building shadow imaging, calculate the value of the satellite altitude angle α.

[0034] Step 2-5, calculate the ratio between the building height and the length of the building shadow The distance and area measurement tools in Google Earth were used to measure the length of building shadows within the target area of ​​the image where a 3D map was to be constructed, and the actual height of the building was estimated based on H = δ·L.

[0035] Furthermore, the step 2-2 includes the following process:

[0036] First, the image is grayscaled to generate a grayscale image; then the grayscale image is histogram equalized to enhance the original image; then the Canny edge detection operator is used to identify and detect building shadows in the preprocessed image.

[0037] Furthermore, in step 4, the size of the search area radius r is determined by the following formula:

[0038]

[0039] Among them, GDOP is the geometric dilution of precision.

[0040] Furthermore, in step 7, the satellite S is predicted i At candidate position Cj The visibility at is expressed as follows:

[0041]

[0042] in, Satellite S i Altitude angle, E B is the jth candidate position C j Azimuth in the sky occlusion map The building height angle at 1 indicates satellite S i At candidate position C j Can be seen everywhere, 0 means satellite S i At candidate position C j It is not visible anywhere.

[0043] Furthermore, the step 8 specifically includes: setting the minimum threshold SNR of the LOS signal LOS and the highest threshold SNR of NLOS signals NLOS , judge when the received GNSS signal signal-to-noise ratio is greater than SNR NLOS When the received GNSS signal noise ratio is less than SNR LOS The satellite is not visible.

[0044] Furthermore, the step 9 specifically includes:

[0045] At each candidate location, when step 7 shows that the satellite is visible at the candidate location, the following judgment and scoring are performed:

[0046] When a GNSS signal is received, it is scored according to the signal-to-noise ratio: when the signal-to-noise ratio is greater than or equal to SNR NLOS When the received GNSS signal noise ratio is less than SNR NLOS When the signal-to-noise ratio is greater than SNR NLOS and is less than SNR NLOS When , the score is 0.5;

[0047] When no GNSS signal is received, the score is 0;

[0048] At each candidate location, if the satellite is not visible at the candidate location as determined in step 7, the following judgment and scoring are performed:

[0049] When a GNSS signal is received, it is scored according to the signal-to-noise ratio: when the signal-to-noise ratio is greater than or equal to SNR NLOS When the received GNSS signal noise ratio is less than SNR NLOS When the signal-to-noise ratio is greater than SNRNLOS and is less than SNR NLOS When , the score is 0.5;

[0050] When no GNSS signal is received, the score is 1.

[0051] Furthermore, the step 11 specifically includes the following process:

[0052] Select state quantity X=[x c ,v c ,x a ,v a ], where x c , x a are the position components of the user’s cross-street direction and along-street direction, v c and v a are the velocity components of the user in the direction of crossing the street and along the street, respectively, and the system continuous state equation is established:

[0053] x(t)=Fx(t)+ω(t)

[0054] Where x(t) is the state vector, t is the time variable, F is the state transfer matrix, and ω(t) is the process noise;

[0055] Assuming that the data sampling period is T, the continuous system state equation is discretized to obtain the discrete state equation of the system:

[0056] X k =Φ k-1 X k-1 +W k-1

[0057] Where W k-1 is a discrete noise sequence, Φ k-1 is the state transfer matrix in the discrete time domain, which is derived from the continuous time domain F, and we get:

[0058]

[0059] In the formula, the state noise satisfies the state noise covariance Here we assume that the processing noise of each state variable is independent of each other and is white noise, and its respective power spectrum density is diag(σ p ,σ v ,σ p ,σ v ), where σ p is the position covariance, that is, the noise power spectrum density of the position; σ v is the noise power spectrum density of the velocity; the discretized Q k-1 for:

[0060]

[0061] wherein Q p is the noise covariance of position and velocity,

[0062] P is the shadow matching positioning solution obtained in step 10 SM is the cross-street direction position component is the along-street direction position component and the velocity v output by the GNSS as observation, i.e. The observation equation between the continuous observation and the state variable is as follows:

[0063]

[0064] wherein ε c and ε a are the observation noises of the cross-street and along-street positions of the shadow matching positioning output respectively; δ v is the velocity observation noise output by the GNSS receiver;

[0065] The recursive filtering equation of the system is obtained according to the linearization of the observation equation, the EKF recursive equation and the established system state equation, and finally the optimal estimation value of each state quantity in the state vector is calculated; and then the fused position solution P is calculated according to the cross-street direction position component and the along-street direction position component result (E result , N result ) calculated by the EKF, wherein N result is the east direction coordinate of the fused position solution, and N result is the north direction coordinate of the fused position solution.

[0066] The application further provides a 3D map-based low-cost positioning system in an urban canyon environment, comprising a computer program which implements the steps in the 3D map-based low-cost positioning method in an urban canyon environment.

[0067] Compared with the prior art, the application has the following advantages and beneficial effects:

[0068] (1) The present invention utilizes the free high-precision two-dimensional image data provided by Google Earth. Using shadow altimetry, the height information of high-rise buildings in urban canyon areas is extracted relatively simply and conveniently. This information is then integrated with the building outline information obtained through OSM to construct a 3D map suitable for shadow matching positioning. This effectively addresses the limitations imposed by the limited popularity of high-precision 3D maps on the application of shadow matching positioning. Furthermore, this method offers a low cost for establishing 3D maps. The present invention only requires establishing and storing maps of areas with weak GNSS signals, and the maps only require extracting building boundary points and corresponding height information. High-precision 3D maps are not required, making the operation simple and effectively avoiding the storage and computing pressure placed on terminal devices by storing and processing high-precision 3D maps.

[0069] (2) The present invention introduces GNSS velocity information on the basis of traditional shadow matching positioning, and fuses the shadow matching positioning solution with GNSS velocity information through extended Kalman filtering to obtain the optimal position solution, which can effectively improve the problem of poor positioning accuracy of traditional shadow matching positioning in the cross-street direction.

[0070] (3) The present invention does not require the introduction of other sensors and has no additional hardware costs. In addition, the cost of building and storing maps is also low. Figure 1 Once built, it can be reused continuously without major changes to the building entity. Therefore, the present invention has the advantage of low cost and is suitable for low-cost terminal devices, making it able to effectively improve the positioning accuracy in urban canyon environments based on existing hardware equipment. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 This is a schematic flow chart of the low-cost positioning method for an urban canyon environment based on a 3D map provided by the present invention. DETAILED DESCRIPTION

[0072] The technical solutions provided by the present invention will be described in detail below with reference to specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0073] The low-cost positioning method for urban canyon environments based on 3D maps provided by the present invention utilizes OSM and Google Earth to establish a 3D map of the target area. On this basis, the building information provided by the 3D map is used to perform shadow matching positioning calculations to assist GNSS positioning. The traditional shadow matching positioning method is improved to enhance the positioning accuracy of the shadow matching itself. At the same time, GNSS information is introduced for fusion to improve the overall positioning accuracy.

[0074] The method of the present invention includes an offline map preparation phase and an online positioning phase. The offline map preparation phase aims to establish a 3D map of the urban canyon area for shadow matching. Based on the two-dimensional map, the building height is calculated by using the altitude and azimuth angles of the sun and satellite and the geometric relationship between the image shadow and the object height, thereby producing a local 3D map containing building height information. Specifically, the method includes steps 1 to 3:

[0075] Step 1: Extract building outline data:

[0076] In OpenStreetMap, select the target area for 3D mapping by entering a latitude and longitude range or manually selecting a rectangle, and save the selected area as an .osm file. In ArcMap, the map editing software, use the OpenStreetMap Toolbox to load the saved .osm file and extract the 2D geographic data, including road information and building outlines.

[0077] Step 2, extracting building height information data, includes the following sub-steps:

[0078] 2-1, Acquisition of Google Earth 2D images:

[0079] Run Google Earth, select the shooting time of the selected image in Google Earth's historical images, select the spatial area range corresponding to step 1, and download the image of the selected time and area.

[0080] 2-2, Image processing of Google Earth 2D images:

[0081] Perform image preprocessing on the Google Earth image downloaded in 2-1. First, convert the original image to grayscale to generate a grayscale image. Next, perform histogram equalization on the generated grayscale image to enhance the original image.

[0082] Next, the Canny edge detection operator is used to identify and detect building shadows in the preprocessed image. This facilitates the subsequent estimation of building height and actual shadow length based on the identification and detection of building shadows.

[0083] 2-3, using shadow altimetry to estimate building height, includes the following sub-steps:

[0084] 2-3-1, calculate the solar altitude angle:

[0085] The calculation of the sun elevation angle β is done using the National Oceanic and Atmospheric Administration (NOAA) online tool for calculating the sun elevation angle. Since the two-dimensional images in Google Earth are usually taken in the morning, the shooting time of the image obtained in step 2-1 can be set to 10:00 am. Assuming that the center point of the image area in step 2-1 is P, the sun elevation angle at 10:00 am on the image shooting date is calculated according to the latitude and longitude data of point P and the average elevation of the point.

[0086] 2-3-2, Calculation of satellite elevation angle:

[0087] In the image downloaded in step 2-1, a typical high-rise building Building X with known height data is selected as a reference, and the height of Building X is known to be H. Using the distance and area measurement tool in Google Earth, the length of the shadow of Building X imaged on the satellite image is measured as L. According to the geometric relationship between the satellite, the sun, the building, and the shadow, it can be known that:

[0088] H = L tan a tan β / (tan β - tan a) (1)

[0089] Then the value of the satellite elevation angle a can be obtained:

[0090]

[0091] 2-3-3, High-rise building height estimation:

[0092] In the case where the satellite elevation angle and the sun elevation angle have been calculated in steps 2-3-1 and 2-3-2, the ratio between the building height and the building shadow length is calculated as: Next, the distance and area measurement tool in Google Earth is used again to measure the length of the building shadow in the area of the image where the 3D map needs to be established, and the actual height of the building is estimated according to H = δ L.

[0093] Here are two alternative methods for estimating the height of a high-rise building:

[0094] Method 1: Use the Shadow to Height.sav extension in the envi software. This tool calculates the height of features based on their shadows. Enter the solar altitude and azimuth information in Shadow to Height.sav. The tool automatically calculates the vertical height of the specified feature based on the length of the shadow drawn. This makes it suitable for estimating the height of high-rise buildings based on shadows. In this method, the solar azimuth is measured using the Ruler tool in Google Earth Pro. Use the Ruler tool to add a line from a point in the shadow to the corresponding point on the building. The azimuth displayed by the Ruler tool is the azimuth of the sun at the time the image was captured.

[0095] Method 2: Using Google SketchUp, Google's 3D design software, first create a 3D model of the target building in SketchUp. Then, set the sun's altitude so that the created 3D building model casts a shadow on the ground. Then, adjust the sun's azimuth so that the direction of the 3D model's shadow aligns with the sun's azimuth of the Google Earth 2D image downloaded in 2-1. Next, adjust the 3D model's height so that the shadow cast by the 3D model coincides with the building in the 2D image. At this point, the height of the 3D model can be considered the building's true height. To align the direction of the 3D model's shadow with the sun's azimuth of the Google Earth 2D image downloaded in 2-1, import the 2D image into SketchUp as a basemap. Create a 3D model at the location of the building whose height you want to calculate. Adjust the time in SketchUp so that the shadow cast by the 3D building model in SketchUp coincides with the shadow in the 2D image. This indicates that the sun's azimuth at this time is consistent with the sun's azimuth when the 2D image was captured. Then, use the stretch tool to adjust the height of the 3D building model so that the shadow cast by the software completely overlaps with the shadow in the 2D image. At this point, use a tape measure to measure the height of the 3D model, which is the actual building height.

[0096] Step 3: Fusion of 2D building outlines and building height information:

[0097] Import the building height information extracted in step 2 into the two-dimensional geographic data extracted in step 1 to generate a 3D map containing building height information suitable for shadow matching positioning.

[0098] The online positioning phase uses the 3D map prepared in the offline phase, an improved shadow matching positioning algorithm, and GNSS velocity information to locate the terminal device. Specifically, the following steps are included: 4-11:

[0099] Step 4: Determine the search area:

[0100] GNSS raw signal data is received using GNSS receivers, mobile phones, and other terminal devices. Pseudorange information and broadcast ephemeris are used to perform pseudorange single-point positioning. The coordinates of the initial position P in the Earth-centered Earth-fixed coordinate system are obtained and converted to Gaussian projection coordinates. A circular search area with P as the center and r as the radius is generated. The size of the search area radius r is determined by the quality of the satellite data, with the geometric dilution of precision (GDOP) as the evaluation indicator. The size of r is determined by formula (3):

[0101]

[0102] Step 5: Generate candidate location grid:

[0103] Within the circular search area determined in step 4, candidate locations are set at intervals of 1 m to create a candidate location grid.

[0104] Step 6: Generate a sky occlusion map for each candidate location:

[0105] At each candidate location, the building height information is used to generate a sky occlusion map. The sky occlusion map contains the elevation angle information of the corresponding building boundary at each angle within the 360° azimuth range.

[0106] Step 7, predict satellite visibility:

[0107] Get satellite S through broadcast ephemeris i Azimuth and altitude angle Compare the jth candidate position C j Azimuth in the sky occlusion map Building height angle E B and satellite elevation angle The relationship between the predicted satellite S i At candidate position C j Visibility of the location, judgment method:

[0108]

[0109] Prediction satellite S i At candidate position C j The visibility at is expressed as:

[0110]

[0111] Step 8: Observe satellite visibility

[0112] Based on empirical values, the SNR threshold for LOS signals is set to 35dB, and the threshold for NLOS signals is set to 25dB. By comparing the SNR of the received GNSS signal with the set threshold, the type of satellite signal is determined, and the visibility of the observed satellite is determined. The specific judgment method is shown in the following table:

[0113] Signal-to-noise ratio (dB-Hz) Signal Type Satellite visibility >35 LOS visible <25 NLOS Invisible

[0114] In addition, a carrier-to-noise ratio threshold can also be set to determine the type of satellite signal and the visibility of the observed satellite by comparing the carrier-to-noise ratio of the received GNSS signal with the set threshold.

[0115] Step 9, matching scoring:

[0116] At each candidate location, score the match between the predicted satellite visibility in step 7 and the observed satellite visibility in step 8 according to the following scoring template:

[0117]

[0118] Indicates visible, Indicates invisible.

[0119] When the carrier-to-noise ratio threshold is used to determine the visibility of the observed satellite in step 8, this step should be combined with step 7 to perform a comprehensive scoring judgment based on the carrier-to-noise ratio threshold.

[0120] Step 10: Select the candidate position with the highest score to calculate the shadow matching positioning solution:

[0121] After step 9, the total score of all satellites at each candidate position is calculated. The higher the score, the better the match between the predicted satellite visibility and the observed satellite visibility at the candidate position, and the closer the position is to the actual location of the user receiving device. Since different candidate positions may have the same score, the candidate position with the highest score is selected, and the average of the easting and northing coordinates of these candidate positions in the Gaussian plane rectangular coordinate system is taken as the shadow matching positioning solution P SM Coordinates:

[0122]

[0123] where N i and E i are the northing and easting coordinates of the i-th candidate location with the highest score, and k is the number of candidate locations with the highest score. N SM and E SM That is P SM The northing and easting coordinates of .

[0124] Step 11: Fuse the shadow matching positioning solution with the GNSS velocity information:

[0125] Based on the shadow matching positioning solution obtained in step 10, GNSS velocity information is introduced. Through EKF design, the state equation and observation equation of the system are established. The shadow matching positioning solution and GNSS velocity information are fused to obtain the optimal position estimate:

[0126] Select state quantity X=[x c ,v c ,x a ,v a ], where x c , x a are the position components of the user’s cross-street direction and along-street direction, v c and v a are the velocity components of the user in the direction of crossing the street and along the street, respectively, and the system continuous state equation is established:

[0127] x(t)=Fx(t)+ω(t) (6)

[0128] Where x(t) is the state vector, t is the time variable, F is the state transfer matrix, and ω(t) is the process noise.

[0129] Assuming that the data sampling period is T, the continuous system state equation is discretized to obtain the discrete state equation of the system:

[0130] X k =Φ k-1 X k-1 +W k-1 (7)

[0131] Where W k-1 is a discrete noise sequence, Φ k-1 is the state transfer matrix in the discrete time domain, which is derived from the continuous time domain F, and we get:

[0132]

[0133] In the formula, the state noise satisfies the state noise covariance Here we assume that the processing noise of each state variable is independent of each other and is white noise, and its respective power spectrum density is diag(σ p ,σ v ,σ p ,σ v ), where σ p is the position covariance, that is, the noise power spectrum density of the position; σ v is the speed noise power spectrum density. Discretized Q k-1 for:

[0134]

[0135] Where: Q p is the noise covariance of position and velocity. The noise power spectrum density of position σ p =0, the speed noise power spectrum density σ v =10 -4 .

[0136] The single epoch shadow matching algorithm positioning result P obtained in step 10 SM Position component in the cross-street direction and the street direction position component And the velocity v output by GNSS, as the observation quantity, that is The observation equation between continuous observations and state variables is:

[0137]

[0138] Where: ε c and ε a are the observation noises of the cross-street and along-street positions output by shadow matching positioning; ε v is the velocity observation noise output by the GNSS receiver, which can be approximated as Gaussian white noise.

[0139] The observation equation is linearized, and the recursive filter equation of the system can be obtained according to the EKF recursive equation and the established system state equation, and the best estimated value of each state quantity in the state vector is finally calculated. Then, according to the cross-street position component calculated by EKF and the street-direction position component Calculate the coordinates P of the fused position solution in the Gaussian plane rectangular coordinate system result (E result ,N result ), where N result is the east coordinate of the fused position solution, N result is the north coordinate of the fused position solution.

[0140] The method of the present invention can be implemented using software. Therefore, the present invention also provides a low-cost positioning system for an urban canyon environment based on a 3D map, including a computer program that implements the steps of the low-cost positioning method for an urban canyon environment based on a 3D map.

[0141] It should be noted that the above content only illustrates the technical idea of the present application, and cannot limit the protection scope of the present application. For ordinary skilled in the art, without departing from the principle of the present application, a number of improvements and refinements can be made, which fall within the protection scope of the claims of the present application.

Claims

1. A low-cost positioning method for urban canyon environments based on 3D maps, characterized by: The steps include: Step 1: Extract building outline data: Use map software to obtain the point, line, and surface two-dimensional geographic data of the target area for which a 3D map needs to be created; Step 2: Extract building height information data: Obtain a two-dimensional image of the target area in step 1 using Google Earth, perform image processing on the two-dimensional image, measure the length of the building shadow in the processed two-dimensional image, and use shadow height measurement to estimate building height information based on the relationship between building height and building shadow length. This specifically includes the following sub-steps: Step 2-1, running Google Earth, selecting the shooting time of an image in the historical images of Google Earth, selecting the target area range corresponding to step 1, and downloading the 2D image of the selected time and target area; Step 2-2, performing image processing on the downloaded image; Step 2-3, use the NOAA solar altitude angle online calculation tool to complete the calculation of the solar altitude angle β; Steps 2-4: Select a typical high-rise building with a known height of H as a reference. Use the distance and area measurement tools in Google Earth to measure the length of the high-rise building's shadow as L. Based on the principle of high-rise building shadow imaging, calculate the value of the satellite altitude angle α. Step 2-5, calculate the ratio between the building height and the length of the building shadow Use the distance and area measurement tools in Google Earth to measure the length of building shadows within the target area of ​​the image where a 3D map is to be created, and estimate the actual height of the building based on H = δ·L. Step 3: Fusion of 2D building outlines and building height information: Importing the building height information extracted in step 2 into the two-dimensional geographic data extracted in step 1 to generate a 3D map suitable for shadow matching positioning that includes the building height information; Step 4: Determine the search area: Use the terminal device to receive GNSS raw signal data, use pseudorange information and broadcast ephemeris to perform pseudorange single-point positioning, obtain the coordinates of the initial position P in the Earth-centered Earth-fixed coordinate system and convert them into Gaussian projection coordinates, determine the size of the search area radius r based on the quality of the satellite data, and generate a circular search area with P as the center and r as the radius; Step 5: Generate candidate location grid: Within the circular search area determined in step 4, candidate positions are set according to preset intervals to establish a candidate position grid; Step 6: Generate a sky occlusion map for each candidate location: Generate a sky occlusion map at each candidate location using building height information, the sky occlusion map containing the elevation angle information of the building boundary corresponding to each angle within a 360° azimuth range; Step 7, predict satellite visibility: Get satellite S through broadcast ephemeris i Azimuth and altitude angle Compare the jth candidate position C j Azimuth in the sky occlusion map Building height angle E B and satellite elevation angle The relationship between the predicted satellite S i At candidate position C j visibility of the place; Step 8: Observe satellite visibility Set a signal-to-noise ratio (SNR) or carrier-to-noise ratio (CNR) threshold based on empirical values. Determine the type of satellite signal and the visibility of the observed satellite by comparing the SNR or CNR of the received GNSS signal with the set threshold. Step 9, matching scoring: At each candidate location, score the match between the satellite visibility predicted in step 7 and the satellite visibility observed in step 8; Step 10: Select the candidate position with the highest score to calculate the shadow matching positioning solution: Count the total scores of all satellites at each candidate position after step 9, select the candidate position with the highest score, and take the average of the easting and northing coordinates of these candidate positions as the coordinates of the shadow matching positioning solution; Step 11: Fuse the shadow matching positioning solution with the GNSS velocity information: Based on the shadow matching positioning solution obtained in step 10, the GNSS velocity information is introduced, and the state equation and observation equation of the system are established through EKF design. The shadow matching positioning solution and GNSS velocity information are fused to obtain the optimal position estimate.

2. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1 is characterized in that: The point-line-surface two-dimensional geographic data of step 1 is obtained by the following process: selecting a target area for which a 3D map needs to be created in OpenStreetMap, importing the target area file into the map editing software ArcMap to extract the point-line-surface two-dimensional geographic data, wherein the point-line-surface two-dimensional geographic data includes road geographic information and building outlines.

3. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1 is characterized in that: The step 2-2 includes the following process: First, the image is grayscaled to generate a grayscale image; then the grayscale image is histogram equalized to enhance the original image; then the Canny edge detection operator is used to identify and detect building shadows in the preprocessed image.

4. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1, characterized in that: In step 4, the size of the search area radius r is determined by the following formula: Among them, GDOP is the geometric dilution of precision.

5. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1 is characterized in that: In step 7, the satellite S is predicted i At candidate position C j The visibility at is expressed as follows: in, Satellite S i Altitude angle, E B is the jth candidate position C j Azimuth in the sky occlusion map The building height angle at 1 indicates satellite S i At candidate position C j Can be seen everywhere, 0 means satellite S i At candidate position C j It is not visible anywhere.

6. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1, characterized in that: The step 8 specifically includes: setting the minimum threshold SNR of the LOS signal LOS and the highest threshold SNR of NLOS signals NLOS , judge when the received GNSS signal signal-to-noise ratio is greater than SNR NLOS When the received GNSS signal noise ratio is less than SNR LOS The satellite is not visible.

7. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1, characterized in that: The step 9 specifically includes: At each candidate location, when step 7 shows that the satellite is visible at the candidate location, the following judgment and scoring are performed: When a GNSS signal is received, it is scored according to the signal-to-noise ratio: when the signal-to-noise ratio is greater than or equal to SNR NLOS When the received GNSS signal noise ratio is less than SNR NlOS When the signal-to-noise ratio is greater than SNR NLOS and is less than SNR NLOS When , the score is 0.5; When no GNSS signal is received, the score is 0; At each candidate location, if the satellite is not visible at the candidate location as determined in step 7, the following judgment and scoring are performed: When a GNSS signal is received, it is scored according to the signal-to-noise ratio: when the signal-to-noise ratio is greater than or equal to SNR NLOS When the received GNSS signal noise ratio is less than SNR NLOS When the signal-to-noise ratio is greater than SNR NLOS and is less than SNR NLOS When , the score is 0.5; When no GNSS signal is received, the score is 1.

8. The low-cost positioning method for urban canyon environments based on 3D maps according to claim 1, characterized in that: The step 11 specifically includes the following process: Select state quantity X=[x c ,v c ,x a ,v a ], where x c , x a are the position components of the user’s cross-street direction and along-street direction, v c and v a are the velocity components of the user in the direction of crossing the street and along the street, respectively, and the system continuous state equation is established: x(t)=Fx(t)+ω(t) Where x(t) is the state vector, t is the time variable, F is the state transfer matrix, and ω(t) is the process noise; Assuming that the data sampling period is T, the continuous system state equation is discretized to obtain the discrete state equation of the system: X k =Φ k-1 X k-1 +W k-1 Where W k-1 is a discrete noise sequence, Φ k-1 is the state transfer matrix in the discrete time domain, which is derived from the continuous time domain F, and we get: In the formula, the state noise satisfies the state noise covariance Here we assume that the processing noise of each state variable is independent of each other and is white noise, and its respective power spectrum density is diag(σ p ,σ v ,σ p ,σ v ), where σ p is the position covariance, that is, the noise power spectrum density of the position; σ v is the noise power spectrum density of the velocity; the discretized Q k-1 for: Where: Q p is the noise covariance of position and velocity, The shadow matching solution P obtained in step 10 is SM Position component in the cross-street direction and the street-direction position component And the velocity v output by GNSS is used as the observation quantity, that is, The observation equation between continuous observations and state variables is: Where: ε c and ε a are the observation noises of the cross-street and along-street positions output by shadow matching positioning; ε v is the velocity observation noise output by the GNSS receiver; The observation equation is linearized, and the recursive filter equation of the system is obtained according to the EKF recursive equation and the established system state equation. Finally, the best estimated value of each state quantity in the state vector is calculated; then the cross-street position component calculated by EKF is obtained. and the street direction position component Calculate the fused position solution P result (E result ,N result ), where N result is the east coordinate of the fused position solution, N result is the north coordinate of the fused position solution.

9. A low-cost positioning system for urban canyon environments based on 3D maps, comprising a computer program, characterized in that: The computer program implements the steps of the low-cost positioning method for urban canyon environment based on 3D map as described in any one of claims 1 to 8.

Images