Urban canyon environment low-cost positioning method and system based on 3D map

Abstract

The invention provides an urban canyon environment low-cost positioning method and system based on a 3D map. According to the method, high-precision two-dimensional image data provided by Google Earth are utilized, height information of high-rise buildings in an urban canyon area is extracted through a shadow height measurement method, the height information is fused with building contour information obtained through OSM, and a 3D map suitable for shadow matching positioning is constructed; meanwhile, a shadow matching positioning algorithm is improved, shadow matching positioning calculation is carried out by utilizing building information provided by a 3D map, and the shadow matching positioning calculation is fused with GNSS information. According to the invention, the limitation of low popularization degree of a high-precision 3D map on shadow matching positioning application can be effectively overcome, the problem of poor positioning precision of traditional shadow matching positioning in a cross-street direction is improved, and extra hardware cost is avoided; the method is simple in operation, can effectively avoid the storage and calculation pressure brought to terminal equipment by storing and processing a high-precision 3D map, is suitable for low-cost terminal equipment, and can effectively improve the positioning precision in an urban canyon environment.

Description

Technical Field

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

[0002] GNSS (Global Navigation Satellite System) positioning has become an indispensable part of people's daily life and industrial applications. The most common applications are in vehicle navigation systems, smartphone positioning services, map software, etc. In open areas, GNSS pseudo-range single-point positioning can usually achieve meter-level accuracy, which can meet general positioning needs. However, in urban canyon environments with high-rise buildings on both sides, the number of visible satellites drops sharply, and GNSS signals are easily blocked and reflected by high-rise buildings, resulting in signal attenuation, multipath interference, and non-line-of-sight (NLOS) signal reception. The resulting pseudo-range single-point positioning error can reach tens of meters or even no positioning. The introduction of additional sensors such as lidar, cameras, and IMUs can improve positioning accuracy to a certain extent through multi-sensor fusion technology, but there are limitations in practical applications that are restricted by cost, volume, weight, etc. Thanks to the development of 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 errors in urban environments, and provides a new idea for improving GNSS positioning accuracy in urban canyon environments. The traditional shadow matching positioning algorithm uses 3D city map information, satellite ephemeris data and received satellite signal data to match and narrow down the possible area of ​​the receiving device, thereby realizing the positioning of the receiving device. It can effectively improve the positioning accuracy of the street crossing direction in the urban canyon environment, and has the advantages of low cost and small amount of calculation, but the accuracy is poor in the street direction. In addition, the popularity of urban 3D maps with building height information is not high, and the method of using laser radar scanning, total station measurement and other methods to establish 3D maps is costly and complex, which is not conducive to the application of shadow matching algorithm as an auxiliary positioning method in low-cost receiving devices. Shadow matching positioning only needs to use the boundary information of the building, but the existing high-precision 3D city model often contains a lot of additional detailed information, which increases the complexity of map data processing. At the same time, storing a large amount of high-precision model data on low-cost terminal devices will bring additional pressure, affecting the operating efficiency and resource utilization of the equipment. Summary of the invention

[0003] In view of the problems existing in the prior art, the present invention provides a low-cost positioning method and system for an urban canyon environment based on a 3D map. A 3D map of the target area is established through OpenStreetMap (OSM, an open source map project) and Google Earth, and shadow matching positioning calculation is performed using the building information provided by the 3D map to assist GNSS positioning. At the same time, an improved method is proposed for the traditional shadow matching positioning algorithm to improve the positioning accuracy of the shadow matching itself, and to fuse it with the GNSS information to improve the comprehensive positioning accuracy.

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

[0005] The low-cost positioning method for urban canyon environment based on 3D map includes the following steps:

[0006] Step 1, extract building outline data:

[0007] 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 established;

[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, use the shadow height measurement method, and estimate the building height information based on the relationship between the building height and the building shadow length;

[0010] Step 3: Merge the 2D building outline with the 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 containing the building height information;

[0012] Step 4, determine the search area:

[0013] Use the terminal device to receive the GNSS raw signal data, use the 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 according to the satellite data quality, and generate a circular search area with P as the center and r as the radius;

[0014] Step 5, generate candidate location grid:

[0015] In 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 sky occlusion maps for each candidate location:

[0017] Generate a sky occlusion map at each candidate location using building height information, wherein the sky occlusion map includes height 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 the altitude angle , compare the jth candidate position C j Azimuth in the sky occlusion map Building height angle E B and satellite altitude 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] The signal-to-noise ratio or carrier-to-noise ratio threshold is set according to the empirical value, and the type of the satellite signal is judged by comparing the signal-to-noise ratio or carrier-to-noise ratio of the received GNSS signal with the set threshold, and the visibility of the observed satellite is determined;

[0022] Step 9, matching scoring:

[0023] Score the match between the satellite visibility predicted in step 7 and the satellite visibility observed in step 8 at each candidate location;

[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 east and north coordinates of these candidate positions as the coordinates of the shadow matching positioning solution;

[0026] Step 11: Fusion of shadow matching positioning solution and 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 through 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, the step 2 includes the following sub-steps:

[0030] Step 2-1, running Google Earth, selecting the shooting time of the selected image in the historical images of Google Earth, selecting the target area range corresponding to step 1, and downloading the two-dimensional 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] Step 2-4, select a typical high-rise building with a known height data of H as a reference, use the distance and area measurement tool in Google Earth to measure the shadow length of the high-rise building as L, and according to the shadow imaging principle of the high-rise building, obtain 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 Use the distance and area measurement tools in Google Earth to measure the length of the building shadow in the target area of ​​the image where a 3D map needs to be built, and estimate the actual height of the building based on H = δ.L.

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

[0036] First, the image is grayed 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 the 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 The 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 Can't be seen 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 signal NLOS , when the received GNSS signal 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 position, when step 7 shows that the satellite is visible at the candidate position, 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 position, when step 7 shows that the satellite is not visible at the candidate position, 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, respectively, and v c and v a are the speed components of the user in the direction of crossing the street and along the street, respectively, and the continuous state equation of the system is established:

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

[0054] Among them, 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 discrete time domain, which is derived from F in continuous time domain, and we get:

[0058]

[0059] In the formula, the state noise satisfies the state noise covariance It is assumed here 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] Where: Q p is the noise covariance of position and velocity,

[0062] 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 observed quantity, that is, The observation equation between continuous observables and state variables is:

[0063]

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

[0065] 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 position component of the cross-street direction 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.

[0066] The present invention also provides a low-cost positioning system for an urban canyon environment based on a 3D map, including a computer program, which implements the steps in the low-cost positioning method for an urban canyon environment based on a 3D map.

[0067] Compared with the prior art, the present invention 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, and extracts the height information of high-rise buildings in the urban canyon area in a relatively simple and convenient manner through the shadow height measurement method, and integrates it with the building outline information obtained through OSM to construct a 3D map suitable for shadow matching positioning, which can effectively solve the limitation of the low popularity of high-precision 3D maps on the application of shadow matching positioning, and the method has a low cost for establishing 3D maps. The present invention only needs to establish and store maps of weak GNSS signal areas, and the maps only need to extract building boundary points and corresponding height information, without the need for high-precision 3D maps, so the operation is simple, and the storage and computing pressure brought to the terminal device by storing and processing high-precision 3D maps can be effectively avoided.

[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 the 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 does not require additional hardware costs. In addition, the cost of establishing and storing maps is also low. Figure 1 Once constructed, 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, so that it can effectively improve the positioning accuracy in the urban canyon environment based on existing hardware devices. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0072] The technical solution provided by the present invention will be described in detail below in conjunction with specific embodiments. It should be understood that the following specific implementation methods 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 environment based on 3D map provided by the present invention uses OSM and Google Earth to establish a 3D map of the target area, and on this basis uses the building information provided by the 3D map to perform shadow matching positioning calculation to assist GNSS positioning; and improves the traditional shadow matching positioning method to improve the positioning accuracy of the shadow matching itself. At the same time, GNSS information is introduced for fusion to improve the comprehensive positioning accuracy.

[0074] The method of the present invention includes an offline map preparation stage and an online positioning stage, wherein the offline map preparation stage aims to establish a 3D map of an urban canyon area for shadow matching, and on the basis of a two-dimensional map, calculates the building height by using the altitude angle and azimuth angle 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 including steps 1 to 3:

[0075] Step 1, extract building outline data:

[0076] In OpenStreetMap, select the target area for 3D map creation by entering the latitude and longitude range or manually selecting the area, and save the selected area as an osm file. In the map editing software ArcMap, use the OpenStreetMap Toolbox to load the osm file saved in the previous step to extract the point, line and surface 2D geographic data containing road geographic information and building outlines.

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

[0078] 2-1, Google Earth 2D image acquisition:

[0079] Run Google Earth, select the shooting time of the selected image in the historical images of Google Earth, select the spatial area range corresponding to step 1, and download the images 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, grayscale the original image 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 the building shadows in the preprocessed image, so as to further estimate the building height and actual length of the shadow based on the subsequent identification and detection of the building shadow.

[0083] 2-3, using shadow height measurement method to estimate building height, including the following sub-steps:

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

[0085] The calculation of the solar altitude angle β is completed using the National Oceanic and Atmospheric Administration (NOAA) solar altitude angle online calculation tool. 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 preset to 10:00 am. Assuming that the center point of the image area in step 2-1 is P, the solar altitude angle at 10:00 on the day of image shooting is calculated according to the latitude and longitude data of point P and the average altitude of the point.

[0086] 2-3-2, calculate the satellite altitude angle:

[0087] In the image downloaded in step 2-1, select BuildingX, a typical high-rise building with known height data, as a reference, and the height of BuildingX is known to be H. Use the distance and area measurement tools in Google Earth to measure the length of the shadow of BuildingX on the satellite image as L. According to the shadow imaging principle of high-rise buildings and the geometric relationship between satellites, the sun, buildings, and shadows, we can know that:

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

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

[0090]

[0091] 2-3-3, Estimation of the height of high-rise buildings:

[0092] When the satellite altitude angle and the sun altitude angle have been obtained in steps 2-3-1 and 2-3-2, calculate the ratio between the building height and the building shadow length. Next, we use the distance and area measurement tool in Google Earth again to measure the length of the building shadow in the area of ​​the image where the 3D map needs to be built, and estimate the actual height of the building based on H = δ.L.

[0093] Here are two optional methods for estimating the height of high-rise buildings:

[0094] Method 1: Use the Shadow to Height.sav extension tool in the envi software, which can calculate the height of the object using the shadow cast by the object on the ground. Enter the solar altitude angle and solar azimuth information in Shadow to Height.sav. The tool can automatically calculate the vertical height of the specified object according to the projection length drawn by the annotation, so it is suitable for estimating the height of high-rise buildings based on shadows. In this method, the solar azimuth can be measured using the ruler tool of Google Earth Pro. Use the ruler tool to add a line segment from a certain point in the shadow to the corresponding point of the building. The azimuth displayed by the ruler tool is the azimuth of the sun when the image was taken.

[0095] Method 2: Use Google SketchUp, a 3D design software owned by Google. First, use SketchUp software to build a 3D model of the target building. Then set the solar altitude angle so that the built 3D model of the building will produce a shadow on the ground. Then adjust the solar azimuth so that the direction of the shadow of the 3D model is consistent with the solar azimuth of the Google Earth 2D image downloaded in 2-1. Next, adjust the height of the 3D model so that the shadow generated by the 3D model coincides with the building in the 2D image. At this time, the height value of the 3D model can be considered as the real height of the building. The specific method of making the shadow direction of the 3D model consistent with the solar azimuth of the Google Earth 2D image downloaded in 2-1 is as follows: import the 2D image into SketchUp as a base map, build a 3D model at the building whose height needs to be calculated, adjust the time in SketchUp so that the shadow generated by the 3D model of the building built in SketchUp coincides with the shadow on both sides of the 2D image, which means that the azimuth of the sun at this time is consistent with the azimuth of the sun when the 2D image was taken, and then use the stretching tool to adjust the height of the 3D model of the building so that the shadow generated in the software completely coincides with the shadow in the 2D image. At this time, the height of the 3D model measured with the tape measure tool is the actual building height.

[0096] Step 3: Merge the 2D building outline with the building height information:

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

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

[0099] Step 4, determine the search area:

[0100] Use GNSS receivers, mobile phones and other terminal devices 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. Generate a circular search area with P as the center and r as the radius. The size of the search area radius r is determined according to the quality of satellite data, and the geometric precision factor GDOP is used as the evaluation index. The size of r is determined by formula (3):

[0101]

[0102] Step 5, generate candidate location grid:

[0103] In the circular search area determined in step 4, candidate positions are set at intervals of 1 m to establish a candidate position grid.

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

[0105] The sky occlusion map is generated at each candidate location using the building height information. The sky occlusion map contains the height angle information of the building boundary corresponding to each angle within the 360° azimuth range.

[0106] Step 7, predict satellite visibility:

[0107] Get satellite S through broadcast ephemeris i Azimuth And the altitude angle Compare the jth candidate position C j Azimuth in the sky occlusion map Building height angle E B and satellite altitude 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] According to the empirical value, the signal-to-noise ratio threshold of the LOS signal is set to 35dB, and the threshold of the NLOS signal is set to 25dB. By comparing the relationship between the signal-to-noise ratio of the received GNSS signal and the set threshold, the type of the satellite signal is judged 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 may also be set to judge the type of the satellite signal and determine 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, the match between the predicted satellite visibility in step 7 and the observed satellite visibility in step 8 is scored according to the scoring template shown below:

[0117]

[0118] Indicates visible, Indicates invisible.

[0119] When the carrier-to-noise ratio threshold is used in step 8 to determine the visibility of the observed satellite, this step should be combined with step 7 to make 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 counted. The higher the score, the higher 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 position 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 east and north 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, respectively, 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: Fusion of shadow matching positioning solution and GNSS velocity information:

[0125] 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:

[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, respectively, and v c and v a are the speed components of the user in the direction of crossing the street and along the street, respectively, and the continuous state equation of the system is established:

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

[0128] Among them, 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 discrete time domain, which is derived from F in continuous time domain, and we get:

[0132]

[0133] In the formula, the state noise satisfies the state noise covariance It is assumed here 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. The 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 is 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 observed quantity, that is The observation equation between continuous observables and state variables is:

[0137]

[0138] Where: ε c and ε a are the observation noises of the cross-street and along-street locations output by shadow matching positioning; ε v is the velocity observation noise output by the GNSS receiver, which can be approximated as (0, ) is a 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 finally the best estimated value of each state quantity in the state vector is calculated. 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 by 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 in 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 invention and cannot be used to limit the protection scope of the present invention. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications all fall within the protection scope of the claims of the present invention.

Claims

1. A low-cost positioning method for urban canyon environments based on 3D maps, characterized in that: 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 established; Step 2: Extract building height information data: 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, use the shadow height measurement method, and estimate the building height information based on the relationship between the building height and the building shadow length; Step 3: Merge the 2D building outline with the 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 containing the building height information; Step 4, determine the search area: Use the terminal device to receive the GNSS raw signal data, use the 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 according to the satellite data quality, and generate a circular search area with P as the center and r as the radius; Step 5, generate candidate location grid: In 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 sky occlusion maps for each candidate location: Generate a sky occlusion map at each candidate location using building height information, wherein the sky occlusion map includes height 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 the altitude angle Compare the jth candidate position C j Azimuth angle in the sky occlusion map Building height angle E B and satellite altitude angle The relationship between the predicted satellite S i At candidate position C j visibility of the place; Step 8: Observe satellite visibility The signal-to-noise ratio or carrier-to-noise ratio threshold is set according to the empirical value, and the type of the satellite signal is judged by comparing the signal-to-noise ratio or carrier-to-noise ratio of the received GNSS signal with the set threshold, and the visibility of the observed satellite is determined; Step 9, matching scoring: Score the match between the satellite visibility predicted in step 7 and the satellite visibility observed in step 8 at each candidate location; 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 east and north coordinates of these candidate positions as the coordinates of the shadow matching positioning solution; Step 11: Fusion of shadow matching positioning solution and 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 environment based on 3D map 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 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.

3. The low-cost positioning method for urban canyon environment based on 3D map according to claim 1 is characterized in that: The step 2 includes the following sub-steps: Step 2-1, running Google Earth, selecting the shooting time of the selected image in the historical images of Google Earth, selecting the target area range corresponding to step 1, and downloading the two-dimensional 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 β; Step 2-4, select a typical high-rise building with a known height data of H as a reference, use the distance and area measurement tool in Google Earth to measure the shadow length of the high-rise building as L, and according to the shadow imaging principle of the high-rise building, obtain 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 the building shadow in the target area of ​​the image where a 3D map needs to be built, and estimate the actual height of the building based on H = δ.L.

4. The low-cost positioning method for urban canyon environment based on 3D map according to claim 3 is characterized in that: The step 2-2 includes the following process: First, the image is grayed 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 the building shadows in the preprocessed image.

5. The low-cost positioning method for urban canyon environment based on 3D map according to claim 1 is 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.

6. The low-cost positioning method for urban canyon environment based on 3D map 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 The altitude angle, E B is the jth candidate position C j Azimuth angle in the sky occlusion map The building height angle at 1 means satellite S i At candidate position C j Can be seen everywhere, 0 means satellite S i At candidate position C j Can't be seen anywhere.

7. The low-cost positioning method for urban canyon environment based on 3D map according to claim 1 is 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 signal NLOS , when the received GNSS signal noise ratio is greater than SNR NLOS When the received GNSS signal noise ratio is less than SNR LOS The satellite is not visible.

8. The low-cost positioning method for urban canyon environment based on 3D map according to claim 1 is characterized in that: The step 9 specifically includes: At each candidate position, when step 7 shows that the satellite is visible at the candidate position, 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 position, when step 7 shows that the satellite is not visible at the candidate position, 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.

9. The low-cost positioning method for urban canyon environment based on 3D map 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, respectively, and v c and v a are the speed components of the user in the direction of crossing the street and along the street, respectively, and the continuous state equation of the system is established: x(t)=Fx(t)+ω(t) Among them, 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 discrete time domain, which is derived from F in continuous time domain, and we get: In the formula, the state noise satisfies the state noise covariance It is assumed here 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 observed quantity, that is, The observation equation between continuous observables and state variables is: Where: ε c and ε a are the observation noises of the cross-street and along-street locations 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 position component of the cross-street direction 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.

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

Images