3D Grid-based Global Positioning System Satellite Signal Shadowing Range Modeling in Urban Area

In urban environments, the location accuracy of autonomous vehicles and unmanned aerial vehicles (UAVs) using navigation systems may be degraded because the global navigation satellite system (GNSS) positioning accuracy is degraded owing to satellite visibility or multipath signals. Therefore, to ensure the safety of autonomous navigation systems, we need technologies that can recognize or predict changes in the GNSS observation environment when vehicles enter areas with severe reception conditions. In this study, we developed an algorithm determining the GPS signal shadowing range that can be blocked by nearby obstacles at a specific grid point. To evaluate the algorithm, we collected GPS signals in the test area and analyzed the signal characteristics by overlapping the satellite trajectory and satellite signal shadowing area modeling results. As a result of the experiment, GPS signals were received in the satellite invisible area owing to the multipath effect, but the data continuity was reduced in the corresponding satellites, and L2 signals were disconnected or the amplitude of the vibration of MP1 was significantly larger than that of the visible satellites. Therefore, the SEM generation algorithm is expected to improve the positioning accuracy of autonomous navigation systems by removing multipath signals and simulating the observation environment of the GNSS receiver.

In urban environments, the location accuracy of autonomous vehicles and unmanned aerial vehicles (UAVs) using navigation systems may be degraded because the global navigation satellite system (GNSS) positioning accuracy is degraded owing to satellite visibility or multipath signals. Therefore, to ensure the safety of autonomous navigation systems, we need technologies that can recognize or predict changes in the GNSS observation environment when vehicles enter areas with severe reception conditions. In this study, we developed an algorithm determining the GPS signal shadowing range that can be blocked by nearby obstacles at a specific grid point. To evaluate the algorithm, we collected GPS signals in the test area and analyzed the signal characteristics by overlapping the satellite trajectory and satellite signal shadowing area modeling results. As a result of the experiment, GPS signals were received in the satellite invisible area owing to the multipath effect, but the data continuity was reduced in the corresponding satellites, and L2 signals were disconnected or the amplitude of the vibration of MP1 was significantly larger than that of the visible satellites. Therefore, the SEM generation algorithm is expected to improve the positioning accuracy of autonomous navigation systems by removing multipath signals and simulating the observation environment of the GNSS receiver.

Introduction
Autonomous driving systems are composed of body control technologies that reflect the decisions made by artificial intelligence (AI) based on accurately recognizing the driving environment and situation. Precise kinematic positioning technology is a key technology of autonomous driving systems because the location of the vehicle needs to be accurately expressed on high-definition (HD) maps to accurately recognize the driving environment. The sensors used for positioning vehicles include the global navigation satellite system (GNSS) and inertial navigation system (INS), and autonomous vehicles use GNSS/INS positioning systems, which combine these two systems to complement each other's strengths and weaknesses. RADAR, LiDAR, optical sensors, and ultrasonic sensors are also applied to improve safety by supplementing the positioning and situational awareness performance of vehicles.
Similarly to autonomous vehicles, unmanned aerial vehicles (UAVs) also use autonomous navigation systems. The latest UAV technologies focus on efficiency and economic feasibility to achieve the desired purpose with less manpower and resources than conventional technologies, and the technology trend is aimed at autonomy or unmanned applications by introducing autonomous navigation technologies as in the case of autonomous vehicles. UAVs are used at lower altitudes than aircraft, so they must avoid collisions with ground structures. Although obstacle avoidance technologies using ultrasonic or optical sensors are being developed, (1,2) positioning accuracy is also important in the application of UAVs and autonomous vehicles in urban areas with many high-rise buildings. (3) Radišić et al. pointed out that the environmental characteristics of urban areas may affect the surveillance and navigation of UAVs. (4) Zimmermann et al. proposed a method to ensure the accuracy of UAV mapping results by improving the GNSS positioning accuracy of UAVs by applying a 3D model-based shadow matching technology. (5) Studies were also carried out on how to use 3D spatial information to support the situational awareness of UAVs, such as the HD maps used in autonomous vehicles. (6,7) These recent studies imply the need to secure safety for autonomous navigation systems such as autonomous vehicles and UAVs.
Conventional techniques for reducing the GNSS multipath effect include antenna-based mitigation and improved receiver technologies, which use antennas such as choke rings and ground planes, and signal and data processing methods. (8) Recently, the spread of drones and the development of drone-based mapping technologies have made it easier to build highresolution 3D terrain data than in the past, which has led to the development of technologies that apply 3D spatial data in areas such as autonomous vehicles and ground survey planning. (9,10) In Korea, research is being performed on a 3D grid system as a data frame for implementing digital twin and UAV traffic control that can link various data of cities in cyberspace. (11) In addition to the sensors in autonomous navigation systems, such grid-based geospatial information could be used as a way for machines to recognize the real world. In addition, grid-type 3D models can be applied to various applications as they make it easier to reduce the data volume and apply physical simulations than vector-based 3D models with complicated geometric relationships.
In autonomous navigation systems, the positioning in outdoor environments is mainly performed through GNSS. Although the absolute location of the user can be determined by using GNSS, positioning becomes impossible in environments where minimal GNSS satellites cannot be observed, such as tunnels, and the positioning accuracy is degraded in areas where obstacles block signals, such as in urban areas. Providing GNSS satellite signal coverage information reflecting ground obstacles to autonomous navigation systems can contribute to improving the positioning accuracy by removing the satellite signals in the non-line-of-sight (NLOS) area, such as multipath signals, to calculate the location of the vehicle. In addition, the GNSS observation environment models can increase the safety of autonomous driving by improving positioning accuracy.
The line-polygon collision test method was previously developed to simulate the GNSS satellite observation environment by analyzing the line-of-sight (LOS) area between observation points and satellites. (12) However, this method requires 3D vector information, which makes it difficult to identify the complex geometric relationships between the locations of users and satellites, and the information on obstacles such as buildings. (13) Manual work was also necessary for 3D modeling to use 3D vector information to analyze GNSS visibility because we need to use 3D building models fused with terrain information or the values of substituting the altitude information for the building vectors extracted from the digital map. On the other hand, grid-based ground object modeling can easily build high-resolution spatial data even in areas without 3D models by semiautomatically constructing a digital surface model (DSM) using mapping technologies based on drones.
In this paper, we proposed a method to model the GNSS observation environment based on DSM, a grid-based 3D model constructed through drone mapping. In general, GNSS receiver firmware or postprocessing software includes a function to set the elevation mask that removes satellite data below a certain elevation to reduce the multipath effect of low-altitude satellites. In this paper, we developed an algorithm that automatically models the range of GNSS satellite signals (azimuth, elevation) that may be blocked by surrounding obstacles at a point inside the DSM data, and defined this model as the surface elevation mask (SEM). When a point in the DSM area is selected as the target point to generate the SEM, the SEM generation algorithm implements a straight line from the target point using the equation of a straight line in the image for every azimuth and determines the maximum elevation angle among the elevation angles calculated through the elevation values of the checkpoints and the target point on each straight line. In Sect. 2 we explain the SEM generation algorithm and, in Sect. 3, we explain the result of GPS observation data analysis in the testbed area. We also analyzed the SEM modeling results and GPS satellite signal observations to evaluate the modeling results and analyze the GPS data characteristics in areas with severe reception conditions.

Development of SEM Generation Algorithm
The SEM generation algorithm is used to calculate the range of azimuth and elevation angles in which GNSS satellite signals can be blocked with obstacles at specific target points based on the DSM. The SEM generation algorithm is also used to calculate the maximum signal shadowed elevation angle in a specific azimuth direction from the target point and compute omnidirectional maximum signal shadowed elevation angles through an iteration procedure. Figure 1 shows a flowchart of the SEM generation algorithm described in Sect. 2. In this paper, the interval of theta (n) was 1°.
To implement the SEM generation algorithm, we need a grid-based data that reflects the height of buildings and surface elevation information from the ground. We can use the DSM generated by processing drone aerial photos or grid-based urban 3D models as the surface elevation information data. In this study, we used spatial information established by drone mapping the premises of the Korea Institute of Civil Engineering and Building Technology (KICT) to develop and evaluate the SEM generation algorithm. Drone mapping for orthophoto and DSM acquisition was performed with the Sony RX100 M2 camera mounted on a quadrotor UAV. The aerial photographs were taken at 75% sidelap and 80% forward overlap, and six ground control points (GCPs) were installed on the ground. A total of 241 aerial photographs taken from the UAV were processed using Pix4D. The accuracy of the orthophoto was less than 3 cm. (14) The ground sampling distance (GSD) of the orthophoto and DSM generated by drone mapping is 5 cm. In addition, the coordinate system consists of a planar rectangular coordinate with the top of the image facing north. Figures 2(a) and 2(b) show an orthophoto and a DSM constructed by drone mapping, respectively. The target area is about 400 m east and 700 m north, and consists of a flat terrain. There are only a few buildings over 20 m in the target area, but there are many buildings concentrated within the area.
In this paper, the process of generating the SEM in the road area is shown as an example to describe the SEM generation algorithm. For this purpose, we need to select a target point to generate the SEM in the road area. Since it is difficult to distinguish roads from other areas by using only the DSM, we generated a road area layer by digitizing the road area in the orthophoto as shown in Fig. 2(c), where the blue and red areas represent the results of digitizing the building and road areas, respectively. In practice, a separate DMS for the building area is not necessary, but we generated the building area layer to describe the areas where signals are blocked by the azimuth and elevation of the GNSS satellite. Finally, the road area DSM was generated by extracting elevation information included in the road area layer from the DSM.
The SEM generation algorithm is configured to calculate the elevation mask for all azimuths at 1° intervals from 0 to 359° azimuth from the target point after a specific target point i is selected in the road area. To do this, when a specific grid coordinate is given as the target point in the DSM, the algorithm has to express a straight line for all azimuths from this point. Figure 3 shows the relationship between the grid coordinate system, the target point (i), the rotation angle (θ), and the azimuth angle (AZ), where the c-axis represents the column of the grid coordinates and the r-axis represents the row of the grid coordinates. When target point i(r i , c i ) is selected in the grid coordinates, a straight line ij starting from i and moving toward θ can be expressed as shown in Fig. 3, where j represents the grid coordinates of the checkpoint on the straight-line ij, and j ranges from target point i to the area of the DSM. ρ is a normal perpendicular to the straight-line ij starting at the origin of the grid coordinates. The azimuth angle (AZ) at target point i is clockwise starting from the N-axis of Fig. 3, and therefore, it is different from the rotation angle θ system starting from the E-axis of Fig. 3 moving counterclockwise. To determine the grid coordinates through which the straight line moving from target point i(r i , c i ) to the rotation angle θ passes, we can use the equation of a straight line based on the polar coordinates as (15) cos cos Here, the grid coordinates of the target point are given as the initial input values in the SEM generation algorithm and θ is given sequentially during the repeated calculation process with values from 0 to 359° at 1° intervals. Therefore, we can determine ρ according to θ.
Once the straight line from the target point i in the θ direction is determined, we can extract the grid coordinates of the checkpoints to calculate the maximum elevation angle by determining the coordinates of the grids on the DSM through which the straight line passes. At this point, the maximum elevation angle is calculated within the range from the target point i to the edge of the given DSM when the straight line is extended in the θ direction. For example, in terms of a straight-line ij, as shown in Fig. 3, the range of calculation may include the top or rightmost grid of the DSM depending on the position of the target point and the slope of the straight line. If the range of calculation is selected from r i where the target point is located to the top range of the image, we can determine the c-axis coordinates for each checkpoint by using Eq. (2). If c j is outside the range of the DSM, the corresponding section is excluded from the grid that needs to be analyzed. By applying this method according to the characteristics of each quadrant with the target point as its origin, we can determine the coordinates of the grids that need to be calculated. (2) Figure 4 shows an example of determining the maximum elevation angle calculation section for a straight line in a specific azimuth from the target point, which shows the results of implementing straight lines in the calculation section from 0° azimuth at 45° intervals. Here, the blue and red areas represent the building and road areas, respectively, and the green straight lines represent the results of implementing straight lines for each azimuth from the target point.
Once the coordinates of the checkpoints are determined, we need to calculate the elevation angle between the target point and each checkpoint. To calculate the elevation angle, we need the horizontal distance in the real world between the target point and the checkpoint and the difference in elevation between the two points. The horizontal distance D from target point i to checkpoint j is calculated as Here, s is the scale factor according to the DSM resolution, which is used to convert the grid coordinate system into the actual distance. The elevation difference Δh between points i and j is determined as Here, h i and h j represent the elevations at points i and j, respectively. Once D and Δh are determined, we can calculate the elevation angle (EL) as Figure 5 shows how the maximum elevation angle is determined for straight lines for each given azimuth. Figure 5(a) shows the result of determining the maximum elevation angle calculation area by implementing a straight line for 60° azimuth from target point i. In terms of comparing the three points a, b, and c, as shown in Fig. 5(b), the maximum elevation is determined at point a where the elevation angle is calculated to be the greatest. We can determine the elevation mask angle at a specific azimuth angle by sequentially calculating the elevations for each target grid and returning the maximum elevation angle calculated so far. The SEM at target point i can be determined by performing this process for all of the azimuths relative to the target point. Figure 6 shows the skyplot of the results of calculating the SEM using points d and e as the target points. The green dots on the map of Fig. 6 indicate the locations of points d and e. Point d is located near the parking lot where most of the buildings are in the north. Point e has roads in the northeast and southeast directions, and obstacles exist at a distance shorter than point d.
In skyplot, the red line represents the SEM and the gray color shows areas where obstacles can block the GNSS satellite signals. The SEM at point d was zero from azimuth angles of 120 to 270° because there were no obstacles in the south. On the other hand, a relatively high SEM was calculated for azimuth angles in the north. The geographical characteristics of the southern direction outside the area where the DSM data is built indicate an open area environment with some obstacles. Therefore, to apply this algorithm to actual services, further studies are necessary to configure the proper horizontal range of grid data (DSM) used to determine the SEM. Owing to the limited range of the DSM used in this study, we can obtain an SEM that reflects the GNSS observation environment of the area when there are many obstacles near the checkpoint as in the case of point e.

Results of Field Test and Data Analysis
To evaluate the SEM generation algorithm, we examined the GPS reception results in the actual test area according to the observation environment. The test was performed at point e of Fig. 6, where we expect many GPS signals to be blocked by obstacles. Figure 7 shows photos of the field test and the surrounding obstacles. As shown in Fig. 7, the test area consists of an environment where buildings exist in the northwest and southwest directions, which may block GPS signals, while roads exist in the northeast and southeast directions, so there will be fewer obstacles in the east than in the west. Figure 8 shows the topographical characteristics of the test area. The green dot in Fig. 8(a) represents the GPS observation point and the white outer square area represents a 50 m radius from the observation point. The elevation information for this area was extracted from the DSM, and Fig. 8(c) shows this information in three dimensions. As shown in Fig. 8(b), the observation point is located in an area where we expect blocked signals because there are buildings to the west. The elevation of the observation point is about 14 m, and the elevation difference between the buildings adjacent to the west is estimated to be more than 10 m. The buildings in the east direction are farther from the observation point than those in the west, but we also expect blocked signals in the east.
The observation test was performed using a Trimble BD970 receiver module and an AV59 antenna for static observation at 1 s epochs for about 20 min. We performed network RTK at 1 s intervals while recording the observation data and calculated the reference coordinates of the test point by averaging all of the positioning results. The raw data acquired by observation was converted to the RINEX file format and preprocessed to analyze the GPS data through TEQC. (16) We used RTKlib, (17) an open-source GNSS data processing program, to extract observation data to analyze the observation results. In terms of the navigation RINEX file required for data processing, we used the broadcast orbit downloaded from the NASA FTP server. (18) Through RTKlib data processing, we acquired the elevation of satellites recorded in the observation data, azimuth information, the data acquisition status by frequency, and MP1 and MP2 data to check the multipath effects of the acquired signals. MP1 and MP2 are linear combinations of pseudorange and phase observations, calculated as (16)  (Color online) Field test photo (center) and photos of obstacles surrounding the test site.
Here, P 1 , P 2 , φ 1 , and φ 2 are the pseudorange and phase observations of L1 and L2, respectively, and α is a constant term calculated as the ratio of the squares of the frequencies f 1 and f 2 ( f 1 2 /f 2 2 ). The two terms MP1 and MP2 allow us to identify the daily RMS values of the L1 and L2 pseudorange multipaths for P code observations, respectively. Figure 9 shows the skyplot of overlapping the satellite azimuth and elevation information obtained by analyzing the SEM and GPS observation data calculated at the observation point. The dots representing the trajectories of the satellites are in yellow, green, red, and blue, which record L1, L1/L2, L1/L5, and L1/L2/L5, respectively. The frequencies provided by satellites differ by block, and Table 1 shows the attributes of the satellites observed in the test. For example, the trajectory of satellite 12 is green because it is a Block IIR-M satellite, which only provides L1/L2 signals. Block IIF satellites added after Block IIR-M satellites through GPS modernization also provide L5 signals. Among the satellites observed in the test, satellites 10, 25, and 32 were Block IIF satellites.
In Fig. 9, the number recorded next to the satellite trajectory represents the PRN number of each GPS satellite, and the point close to the text is the point recorded earlier. For example, we can see that satellite no. 29 was in the SEM area and gradually moved to an open sky environment. In terms of comparing the observation environment by satellite, satellites 12, 25, and 32 are considered to be unaffected by obstacles during the entire observation time. As shown in Fig. 9, some data from satellites located at elevations lower than the SEM have been recorded, but L2 signals were intermittently blocked at an elevation lower than the SEM. In terms of comparing the observation rate of L1 data, satellites 12, 25, and 32 located at elevations higher than the SEM area showed observation rates higher than 99%, followed by satellites 29, 14, 10, and 31. In the case of satellite 31, most of the signals were blocked and could not be recorded or were the results of recording multipath signals. Figure 10 shows the characteristics of GPS satellite signals over time, and the red squares indicate the sections where satellites were at elevations lower than the SEM. Figure 10(a) shows the results of recording observations by each satellite over time. In terms of satellites 12, 25, and 32, which were at elevations higher than the SEM, most of the observations were collected normally. On the other hand, the signals of satellites that entered elevations lower than the SEM were temporarily disconnected or blocked. In particular, L2 signals were blocked more frequently than other signals. Figure 10(b) shows the results of calculating MP1 by each satellite and the elevation change over time. Since both L1 and L2 observations are necessary to calculate MP1, there are gaps in the MP1 calculation results for satellites 10, 14, 29, and 31, owing to blocked signals.  We can see that consecutive MP1 values were calculated for satellites 12, 25, and 32, and satellite 12 showed a large amplitude of MP1 values in the latter half of the observation compared with satellites 25 and 32. This is because the elevation angle of satellite 12 was lower than those of satellites 25 and 32, and the stainless steel constituting the outer wall of the vehicle repair building in the corresponding azimuth caused multipath effects. This phenomenon has also been found in the long-term data analysis of the permanent GPS station near the stainlesssteel structure. (19) The MP1 and MP2 RMS values of satellites 25 and 32, which were at elevations higher than the SEM, were lower than those of other satellites, providing stable observations.

Conclusions and Future Challenges
In this study, we developed an algorithm that can determine the range of GNSS signals that may be blocked by surrounding obstacles at a certain grid point on the ground by using the DSM constructed by drone mapping. As a result of evaluating the algorithm by performing a GPS observation test in an actual target area, we found observations even for GPS satellites within the SEM range, but the data continuity was degraded. The L2 signals were disconnected and the amplitude of the MP1 values was significantly larger than that of the visible satellites, so these signals were considered to be multipath signals. Therefore, in this study, we confirmed that we can model GNSS observation environments almost equivalent to the real world by using a high-resolution DSM. Since obstacle information such as buildings does not change significantly except for special circumstances such as construction sites, the algorithm developed in this study can be easily applied by constructing a DSM by UAV mapping urban areas where the positioning accuracy is degraded to implement services for autonomous navigation systems in the future. In this study, we could not confirm the effects of high-rise buildings outside the DSM area owing to a limited DSM range, so further research is necessary on setting the optimal calculation range to determine the SEM to apply this technology to actual cities. GNSS satellite orbits are flexible but have sidereal periodicity, and we have technologies to predict the orbits. To apply this technology as a method to filter multipath signals, techniques that can accurately match the prediction results of the satellite orbit to the target point are necessary.
Recently, digital twin technologies that connect the real world and the cyber world have been proposed as a solution to implement smart cities. 3D grid-based spatial information is receiving much attention as a means to effectively implement city-level digital twins. This algorithm can be used as a base technology to support digital twin and autonomous navigation systems in the future, as it can be applied to raster data and surface models of 3D grid-based spatial information systems based on voxel and octree.