Analysis of the Possibilities of Using Di ﬀ erent Resolution Digital Elevation Models in the Study of Microrelief on the Example of Terrain Passability

: In this article, we discuss issues concerning the development of detailed passability maps, which are used in the crisis management process and for military purposes. The paper presents the authorial methodology of the automatic generation of these maps with the use of high-resolution digital elevation models (DEMs) acquired from airborne laser scanning (light detection and ranging (LIDAR)) and photogrammetric data obtained from unmanned aerial vehicle (UAV) measurements. The aim of the article is to conduct a detailed comparison of these models in the context of their usage in passability map development. The proposed algorithm of map generation was tested comprehensively in terms of the source of the used spatial data, the resolution, and the types of vehicles moving in terrain. Tests were conducted on areas with a diversiﬁed landform, with typical forms of relief that hinder vehicle movement (blu ﬀ s and streams). Due to the huge amount of data to be processed, the comprehensive analysis of the possibilities of using DEMs in di ﬀ erent conﬁgurations of pixel size was executed. This allowed for decreasing the resolution of the model while maintaining the appropriate accuracy properties of the resulting passability map. The obtained results showed insigniﬁcant disparities between both sources of used DEMs and demonstrated that using the model with the 2.5 m pixel size did not signiﬁcantly degrade the accuracy of the passability maps, which has a huge impact on their generation time.


Introduction
Effective conduction of the action connected with, in the broad sense, crisis management requires appropriate planning. One of the most crucial aspects, which should be considered, is the availability and, consequently, the passability of the terrain, particularly when there is an indispensability of conducting operations on terrains containing no roads. Beyond crisis management, the modelling of the terrain passability for different kinds of vehicles is immeasurably important in the planning of military operations, where the necessity of conducting fast and enemy-surprising maneuvers forces the relocation of military vehicles outside the regular road network.
In accordance with the military standardization documents [1], passability is understood as the possibility of overcoming terrain by vehicles (including specialistic ones, e.g., fire engines and combat vehicles) in all weather conditions, both on roads and cross-country. The process of the passability terrain determination itself, in its easiest form, is based on the division of terrain into three classes: GO, SLOW GO, and NO GO terrain [2]. The development of such maps is, therefore, an example of the terrain classification. We see this most often in the context of the classification of satellite images to extract particular elements of land cover [3], soils [4], and land development [5] classes.

Area of Research
In the passability analysis, three parts of the terrain were used. On two of them, there were microrelief objects that we can visually identify from the orthophotomap or DEM visualization, which are likely to be identified as impassable areas. The third part is a flat terrain, which is a passable terrain from a military point of view, as there are no height differences and other terrain obstacles. This was used to check whether the script works correctly. The test areas of interest were squares measuring 50 m by 50 m, and they are located in the south part of Poland near Tylicz city-49°24′ N, 21°01′ E (Figure 1). The test areas are shown in Figure 2, where the characteristic microrelief shapes are marked. These terrain objects can potentially hinder the movement of vehicles. (1) (2)

Area of Research
In the passability analysis, three parts of the terrain were used. On two of them, there were microrelief objects that we can visually identify from the orthophotomap or DEM visualization, which are likely to be identified as impassable areas. The third part is a flat terrain, which is a passable terrain from a military point of view, as there are no height differences and other terrain obstacles. This was used to check whether the script works correctly. The test areas of interest were squares measuring 50 m by 50 m, and they are located in the south part of Poland near Tylicz city-49°24′ N, 21°01′ E (Figure 1). The test areas are shown in Figure 2, where the characteristic microrelief shapes are marked. These terrain objects can potentially hinder the movement of vehicles. (1) (2)

Geomatics Data Used in the Experimentation
In the conducted research, two sources of data were used. The first was light detection and ranging (LIDAR) data that were obtained from the Head Office of Geodesy and Cartography from National Geodesy and Cartography Resource (denoted later in the article as LIDAR). This is a repository of all analogue and digital materials concerning geodesy and cartography that are used in different fields, such as the national economy, national defence, science, and culture [40]. These data are available for the whole of Poland's territory. Their incontrovertible advantage is the ease in acquisition. We can download a part of LIDAR data for any area located in Poland with use of a web browser. The downloaded data are distributed as text files that contain coordinates (X, Y, Z) of points disposed evenly as a net with spacing of 0.5 m for urban areas (standard II) or 1 m for other areas (standard I). The data were interpolated from the point cloud acquired from the airborne laser scanning (LIDAR). A model with the resolution of 0.5 by 0.5 m was available for all three test areas.
The main principles of LIDAR technology are described in [41,42]. LIDAR data are binary files containing point clouds in LAS (Laser) format, which is designed for the processing of LIDAR point cloud data. The accuracy of the DEM created from this source can reach about 0.03 m (paved surface), 0.1 m (natural surface without vegetation), 0.2 m (surface covered with vegetation), and 0.5 m (very thick vegetation) [43]. According to the official data distributor, the mean squared error of the elevations used for research is within 0.2 m [40].
The second source of used data was the unmanned aerial vehicles photogrammetry data (denoted later as UAV). The use of UAVs allowed us to obtain the low altitude images and also photogrammetry data of the chosen terrain [44]. This data can be used, for example, for generating high resolution DEMs and 3D models of different objects or the point cloud data. The DEMs used in this research were generated from the point cloud data obtained from the aero triangulation process. The data was classified and filtered first. The accuracy of the obtained models was close to the accuracy of the LIDAR data [45] because, from the point cloud, we are able to obtain elevation models with horizontal accuracy of 0.05 m and vertical accuracy of 0.1 m as was shown, e.g., in [46]. Obtaining image data using UAV photogrammetry methods, apart from the obvious advantages (low cost, speed, and accuracy of the obtained data), depends to a large extent on the weather conditions [47,48].
The source data for this research were acquired in February in 2019 year with use of a Trimble UX-5 HP airframe, which was equipped with a Sony a7R camera with a 36.64 mm of focal length. The flying platform was equipped with a dual-frequency GPS receiver recording data with 10 Hz frequency. The camera settings were defined in the manual mode, and the focus of the lens was set to infinity. The flight was executed at a 150 m altitude in the east-west direction, assuming that endlap and sidelap were symmetrical and totalled 75%.
Endlap is the longitudinal overlap (along the direction of flight line) and sidelap is the transversal overlap (see [49,50]). This allowed us to obtain a ground sampling distance (GSD) of 0.02 m. This produced JPEG images, which were processed in Inpho UASMaster software. On the area of interest, there were 10 signalled ground control points (GCP) that were designed and measured and four independent check points (CP). All points were measured with use of real time kinematic (RTK)

Geomatics Data Used in the Experimentation
In the conducted research, two sources of data were used. The first was light detection and ranging (LIDAR) data that were obtained from the Head Office of Geodesy and Cartography from National Geodesy and Cartography Resource (denoted later in the article as LIDAR). This is a repository of all analogue and digital materials concerning geodesy and cartography that are used in different fields, such as the national economy, national defence, science, and culture [40]. These data are available for the whole of Poland's territory. Their incontrovertible advantage is the ease in acquisition. We can download a part of LIDAR data for any area located in Poland with use of a web browser. The downloaded data are distributed as text files that contain coordinates (X, Y, Z) of points disposed evenly as a net with spacing of 0.5 m for urban areas (standard II) or 1 m for other areas (standard I). The data were interpolated from the point cloud acquired from the airborne laser scanning (LIDAR). A model with the resolution of 0.5 by 0.5 m was available for all three test areas.
The main principles of LIDAR technology are described in [41,42]. LIDAR data are binary files containing point clouds in LAS (Laser) format, which is designed for the processing of LIDAR point cloud data. The accuracy of the DEM created from this source can reach about 0.03 m (paved surface), 0.1 m (natural surface without vegetation), 0.2 m (surface covered with vegetation), and 0.5 m (very thick vegetation) [43]. According to the official data distributor, the mean squared error of the elevations used for research is within 0.2 m [40].
The second source of used data was the unmanned aerial vehicles photogrammetry data (denoted later as UAV). The use of UAVs allowed us to obtain the low altitude images and also photogrammetry data of the chosen terrain [44]. This data can be used, for example, for generating high resolution DEMs and 3D models of different objects or the point cloud data. The DEMs used in this research were generated from the point cloud data obtained from the aero triangulation process. The data was classified and filtered first. The accuracy of the obtained models was close to the accuracy of the LIDAR data [45] because, from the point cloud, we are able to obtain elevation models with horizontal accuracy of 0.05 m and vertical accuracy of 0.1 m as was shown, e.g., in [46]. Obtaining image data using UAV photogrammetry methods, apart from the obvious advantages (low cost, speed, and accuracy of the obtained data), depends to a large extent on the weather conditions [47,48].
The source data for this research were acquired in February in 2019 year with use of a Trimble UX-5 HP airframe, which was equipped with a Sony a7R camera with a 36.64 mm of focal length. The flying platform was equipped with a dual-frequency GPS receiver recording data with 10 Hz frequency. The camera settings were defined in the manual mode, and the focus of the lens was set to infinity. The flight was executed at a 150 m altitude in the east-west direction, assuming that endlap and sidelap were symmetrical and totalled 75%.
Endlap is the longitudinal overlap (along the direction of flight line) and sidelap is the transversal overlap (see [49,50]). This allowed us to obtain a ground sampling distance (GSD) of 0.02 m. This produced JPEG images, which were processed in Inpho UASMaster software. On the area of interest, there were 10 signalled ground control points (GCP) that were designed and measured and Remote Sens. 2020, 12, 4146 6 of 32 four independent check points (CP). All points were measured with use of real time kinematic (RTK) technique in a global navigation satellite system (GNSS). The terrain coordinates of the GCPs were determined with mean errors totalling: m x = 0.04 m, m y = 0.04 m, and m z = 0.08 m, whereas in the CPs there were the following mean errors: m x = 0.04 m, m y = 0.05 m, and m z = 0.09 m. The σ 0 error for the whole block was 7.2 µm/1.5 px. As far as the exterior orientation is concerned, their mean standard deviation of translations and rotations are presented in the Table 1. The resulting accuracy of the obtained point cloud was for horizontal coordinates (X, Y) of 0.06-0.07 m and for the vertical coordinate (Z) of 0.1 m. The flight path and the tie points distribution are presented in the Figure A1 in Appendix A.

Preparing and Processing Data
To determine the influence of the pixel size of the used DEMs on the accuracy of the obtained passability maps, analyses were conducted on different resolution DEMs. They were created from two sources of data: LIDAR and the UAV photogrammetry data in four different resolutions (0.5, 1, 2.5, and 5 m). Their visualizations are presented in Figure 3. The aim of the differentiation of pixel sizes was to calculate the processing time of the script, which allowed us to determine the dependence between the running time of the tool and the detail of the model. The increase in the DEM pixel size therefore allowed us to reduce the generation time of the passability map and, consequently, to develop it on a larger area in a shorter time.

Preparing and Processing Data
To determine the influence of the pixel size of the used DEMs on the accuracy of the obtained passability maps, analyses were conducted on different resolution DEMs. They were created from two sources of data: LIDAR and the UAV photogrammetry data in four different resolutions (0.5, 1, 2.5, and 5 m). Their visualizations are presented in Figure 3. The aim of the differentiation of pixel sizes was to calculate the processing time of the script, which allowed us to determine the dependence between the running time of the tool and the detail of the model. The increase in the DEM pixel size therefore allowed us to reduce the generation time of the passability map and, consequently, to develop it on a larger area in a shorter time.   The microrelief shapes are characterized by the small size; thus, to capture them, the research was conducted on a high-resolution DEM. We can assume that the pixel size of the model should not be greater than 5 m, because smaller microrelief objects that could be an obstacle to a vehicle would not be identified. From the LIDAR point cloud, we are able to produce models with a 0.5 m pixel size and vertical accuracy to 0.2 m; therefore, the use of this kind of data is legitimate [40].
To summarize, we acquired two points clouds from different sources, and the following step is to create digital elevation models from the classified points. To do so, we interpolated point data to obtain the pixel value of the DEM. The whole process was conducted in ArcGIS software. In the used tool "LAS Dataset to Raster" there were two types of interpolation: • Binning-The determination of the pixel value relies on examining the points that fall within the defined pixel. The final value can be chosen as the minimum, maximum, or average of the points [51].

•
Triangulation-To determine the final value, the Delaunay triangulation was used. This creates a surface made from a network of triangles and then rasterizes it.
The binning method is mainly used for high-density data, when we are sure that, in each pixel, there will be at least one point [52]. For the digital elevation model, the best setting is to use the minimum value, because the algorithm chooses then the lowest point (bare ground) inside the pixel [53]. However, triangulation is used with low-density data, because the pixel value is determined by the surface of the triangles. To see the difference between these two interpolation methods, a comparison was conducted. The two DEMs used the same parameters; however, different interpolation methods were subtracted. The obtained results are presented in Table 2.  The microrelief shapes are characterized by the small size; thus, to capture them, the research was conducted on a high-resolution DEM. We can assume that the pixel size of the model should not be greater than 5 m, because smaller microrelief objects that could be an obstacle to a vehicle would not be identified. From the LIDAR point cloud, we are able to produce models with a 0.5 m pixel size and vertical accuracy to 0.2 m; therefore, the use of this kind of data is legitimate [40].
To summarize, we acquired two points clouds from different sources, and the following step is to create digital elevation models from the classified points. To do so, we interpolated point data to obtain the pixel value of the DEM. The whole process was conducted in ArcGIS software. In the used tool "LAS Dataset to Raster" there were two types of interpolation: • Binning-The determination of the pixel value relies on examining the points that fall within the defined pixel. The final value can be chosen as the minimum, maximum, or average of the points [51].

•
Triangulation-To determine the final value, the Delaunay triangulation was used. This creates a surface made from a network of triangles and then rasterizes it.
The binning method is mainly used for high-density data, when we are sure that, in each pixel, there will be at least one point [52]. For the digital elevation model, the best setting is to use the minimum value, because the algorithm chooses then the lowest point (bare ground) inside the pixel [53]. However, triangulation is used with low-density data, because the pixel value is determined by the surface of the triangles. To see the difference between these two interpolation methods, a comparison was conducted. The two DEMs used the same parameters; however, different interpolation methods were subtracted. The obtained results are presented in Table 2. As we can see in Table 2, the differences between the two types of interpolations were not high, and, due to the high density of the data, the binning method was chosen.

Terrain Passability Analysis Assumptions
The determination of the terrain passability was conducted for three vehicles used in the Polish Armed Forces: Honker, Star 266, and HMMWV (High Mobility Multipurpose Wheeled Vehicle; colloquially: Humvee). Their appearance is presented in Figure 4. As we can see in Table 2, the differences between the two types of interpolations were not high, and, due to the high density of the data, the binning method was chosen.

Terrain Passability Analysis Assumptions
The determination of the terrain passability was conducted for three vehicles used in the Polish Armed Forces: Honker, Star 266, and HMMWV (High Mobility Multipurpose Wheeled Vehicle; colloquially: Humvee). Their appearance is presented in Figure 4. A Honker is a Polish off-road passenger vehicle, a Star 266 is an off-road truck designed to transport cargo or people on roads and off-roads, whereas a Humvee is a four-wheel drive armoured vehicle that is used in military operations for multiple purposes (from medical to combat). There are certain crucial tactical and technical parameters that influence the cross-country movement of each of the vehicles the most, which will be used in the presented methodology. They are listed in Table  3.  A Honker is a Polish off-road passenger vehicle, a Star 266 is an off-road truck designed to transport cargo or people on roads and off-roads, whereas a Humvee is a four-wheel drive armoured vehicle that is used in military operations for multiple purposes (from medical to combat). There are certain crucial tactical and technical parameters that influence the cross-country movement of each of the vehicles the most, which will be used in the presented methodology. They are listed in Table 3. The passability analysis was divided into two parts. The first part was conducted in QGIS software, where the slope map was created and then the change of slope was calculated. The computed values were compared with the approach and descent angles to identify impassable areas. The impassable areas due to the slope angle are presented in Figure 5. For the analysis, we took the lowest value from both approach and descent angles to include all possible directions of movement.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 32 The passability analysis was divided into two parts. The first part was conducted in QGIS software, where the slope map was created and then the change of slope was calculated. The computed values were compared with the approach and descent angles to identify impassable areas. The impassable areas due to the slope angle are presented in Figure 5. For the analysis, we took the lowest value from both approach and descent angles to include all possible directions of movement. As shown in [36], this method does not detect places where the vehicle can hit the ground with the area of the chassis between the wheels. Thus, to solve this problem, a tool in the Python programming environment was developed. The main working principle of this script was to calculate the chassis plane equation and check if, at a given point, the height of the terrain is greater than the height of the chassis. When such a case occurred, the analyzed point was classified as impassable. The input elevation data was the DEM in ASCII XYZ Grid format. This is a text file that stores the coordinates (X, Y, Z) of regularly spaced points. The size of the grid is dependent on the resolution of the elevation model. The research was conducted iteratively individually for each point of the grid. We assumed that the currently analyzed point of the model is a centroid of the chassis. The script only considered the section of the chassis between wheels what we can see in Figure 6. As shown in [36], this method does not detect places where the vehicle can hit the ground with the area of the chassis between the wheels. Thus, to solve this problem, a tool in the Python programming environment was developed. The main working principle of this script was to calculate the chassis plane equation and check if, at a given point, the height of the terrain is greater than the height of the chassis. When such a case occurred, the analyzed point was classified as impassable. The input elevation data was the DEM in ASCII XYZ Grid format. This is a text file that stores the coordinates (X, Y, Z) of regularly spaced points. The size of the grid is dependent on the resolution of the elevation model. The research was conducted iteratively individually for each point of the grid. We assumed that the currently analyzed point of the model is a centroid of the chassis. The script only considered the section of the chassis between wheels what we can see in Figure 6 Simplistically, the developed algorithm can be characterized as an iterative situating the vehicle model in each point of the elevation model of the analyzed terrain (assuming tangency of its tires with the ground). If, as a result of this positioning, the terrain crosses the chassis, it is assumed that the analyzed part of the terrain is classified as impassable. The research was conducted for different directions of vehicle movement (different angles of vehicle rotation towards the terrain). The developed methodology can be used for different kinds of vehicles. We simply need to know the vehicles' parameters, as listed in Table 3 (then, we are able to change these parameters in the script).

The Principle of Working of the Developed Algorithm
The procedure to determine whether the chassis hits the ground was based on the assumptions of coordinate geometry in three-dimensional space. The movement of a vehicle is considered in four main directions: south-north, east-west, north-east-south-west, and north-west-south-east, and so there are four scripts, which differ in the disposition of the vehicle. It does not matter if the vehicle is moving, for instance, northward or southward because the procedure considers only the wheelbase, track width, and ground clearance (without front and rear overhang, which are typically different) and, thus, may be considered as one case. The working principle of the procedure will be shown on the basis of the south-north direction case. The steps of the process are listed below: 1. The elevation model in ASCII XYZ Grid format is converted into a two-dimensional array with three columns (X, Y, Z) and a number of rows equal to the number of points in the model. The script is, thereby, able to perform the mathematical operations on the data. 2. The procedure from this point iterates through each point of the model (row in the array). The coordinates of selected points are assigned to the centroid of the chassis. On the basis of the wheelbase and track width there are calculated horizontal coordinates of the four wheels. The computation considers the local slope of the terrain because, when the vehicle is pitched, the horizontal dimensions must be reduced. 3. In the next step, the z-coordinates of wheels are calculated. The script chooses four points of the grid that are closest to the wheel and then computes the height with the use of an inverse distance weighted (IDW) interpolation [58]. We can see this step in Figure 7. Simplistically, the developed algorithm can be characterized as an iterative situating the vehicle model in each point of the elevation model of the analyzed terrain (assuming tangency of its tires with the ground). If, as a result of this positioning, the terrain crosses the chassis, it is assumed that the analyzed part of the terrain is classified as impassable. The research was conducted for different directions of vehicle movement (different angles of vehicle rotation towards the terrain). The developed methodology can be used for different kinds of vehicles. We simply need to know the vehicles' parameters, as listed in Table 3 (then, we are able to change these parameters in the script).

The Principle of Working of the Developed Algorithm
The procedure to determine whether the chassis hits the ground was based on the assumptions of coordinate geometry in three-dimensional space. The movement of a vehicle is considered in four main directions: south-north, east-west, north-east-south-west, and north-west-south-east, and so there are four scripts, which differ in the disposition of the vehicle. It does not matter if the vehicle is moving, for instance, northward or southward because the procedure considers only the wheelbase, track width, and ground clearance (without front and rear overhang, which are typically different) and, thus, may be considered as one case. The working principle of the procedure will be shown on the basis of the south-north direction case. The steps of the process are listed below: 1.
The elevation model in ASCII XYZ Grid format is converted into a two-dimensional array with three columns (X, Y, Z) and a number of rows equal to the number of points in the model. The script is, thereby, able to perform the mathematical operations on the data. 2.
The procedure from this point iterates through each point of the model (row in the array). The coordinates of selected points are assigned to the centroid of the chassis. On the basis of the wheelbase and track width there are calculated horizontal coordinates of the four wheels.
The computation considers the local slope of the terrain because, when the vehicle is pitched, the horizontal dimensions must be reduced.

3.
In the next step, the z-coordinates of wheels are calculated. The script chooses four points of the grid that are closest to the wheel and then computes the height with the use of an inverse distance weighted (IDW) interpolation [58]. We can see this step in Figure 7.

4.
Having calculated the coordinates of the wheels, the script computes the equation of the plane containing three tire contact points with the ground. At the beginning, four plane equations are calculated (containing three out of four tangent points in different configurations), and then the procedure searches in which plane position the fourth point of the vehicle (which does not belong to the plane) has the least difference of heights with the theoretical point with the same horizontal coordinates that belongs to the plane. After this test, the best fitted plane is chosen for further analysis. If the difference of heights is greater than the ground clearance, the analyzed point of the elevation model is assigned the value 0, which indicates that this point is impassable. The finding of the plane equation comes to the determination of four parameters of the point normal form of the plane (1): where A, B, C, and D are the parameters. The first three are calculated from the cross product of vectors created between two pairs of analyzed points, The D parameter is computed by the substitution of the coordinates of one analyzed point into formula (2).
where x, y, and z are the coordinates of one of the tangent points. 5.
The aim of the next step is to compute the plane that approximates the chassis plane. To do that, we must calculate the parallel plane to the one that has already been determined. The value of the shift (distance between two planes) is equal to the ground clearance. This situation is shown in Figure 8. 6.
This step extracts the part of the chassis plane that is limited by the wheels. Then, inside the resulting rectangle, the algorithm creates evenly distributed control points. It divides the distances between the wheels (wheelbase and track width) into even segments and places control points between them. As a result, we obtain a grid of points located in the chassis plane, which we can see in Figure 9. The distances between the points depend on the vehicle and especially on their parameters. The algorithm divides the wheelbase and track width values in each vehicle into four even parts. 7.
Finally, the script checks each control point for whether its height exceeds the height of the terrain interpolated from the model. Unless all points meet this requirement (are higher than points on terrain), the analyzed point of the model is assigned the value 0 (impassable). 8.
At the end of the process, the script creates a text file with coordinates of each point and an appended value of the passability (NO GO terrain-0 and GO terrain-1). On the basis of this file, we can distinguish the impassable terrain and calculate the statistics, which are presented further in the later section. Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 32   For the better visualization of the whole process, a block diagram was created. This shows the detailed process that the ASCII GRID XYZ file must undergo to obtain the file with points of the model and the indicator of the passability. The diagram is presented in Figure 10.   For the better visualization of the whole process, a block diagram was created. This shows the detailed process that the ASCII GRID XYZ file must undergo to obtain the file with points of the model and the indicator of the passability. The diagram is presented in Figure 10.   For the better visualization of the whole process, a block diagram was created. This shows the detailed process that the ASCII GRID XYZ file must undergo to obtain the file with points of the model and the indicator of the passability. The diagram is presented in Figure 10. For the better visualization of the whole process, a block diagram was created. This shows the detailed process that the ASCII GRID XYZ file must undergo to obtain the file with points of the model and the indicator of the passability. The diagram is presented in Figure 10. The processing of the vehicle location for each point of the used DEM makes the whole process time consuming. Due to this, to speed up the execution of the script, multiprocessing was set. This indicates that the program evenly divides the whole sum of the elevation model points into several processes, and they are executed simultaneously. Performance tests were conducted on a computer with the following parameters: • Model: ASUS X556UAK • Processor: Intel ® Core™ i5-7200U CPU @ 2.71 GHz • Operating system: Windows 10 ×64 operating system • RAM: 20 GB The obtained execution times of the script on a computer with the above parameters are listed further in the later section. Based on the combination of the analysis results and performance tests, we were able to select the most optimized settings for the whole passability analysis process. This means that we could select the number of processes and the pixel size of the model to obtain the best solution in the fastest way. The results of the tests are described in the next section.

Results
The products of the analysis show the impassable parts of the tested terrains in multiple configurations for each vehicle. The results are computed for each direction of movement (north, east, northeast, and northwest), for two kinds of data (from LIDAR and from UAV) and for each resolution of DEM (0.5, 1, 2.5, and 5 m). Due to the large number of results, some of them are only presented briefly ( Figure 11). The arrow indicates the direction of the vehicle movement. The processing of the vehicle location for each point of the used DEM makes the whole process time consuming. Due to this, to speed up the execution of the script, multiprocessing was set. This indicates that the program evenly divides the whole sum of the elevation model points into several processes, and they are executed simultaneously. Performance tests were conducted on a computer with the following parameters: • Model: ASUS X556UAK • Processor: Intel ® Core™ i5-7200U CPU @ 2.71 GHz • Operating system: Windows 10 ×64 operating system • RAM: 20 GB The obtained execution times of the script on a computer with the above parameters are listed further in the later section. Based on the combination of the analysis results and performance tests, we were able to select the most optimized settings for the whole passability analysis process. This means that we could select the number of processes and the pixel size of the model to obtain the best solution in the fastest way. The results of the tests are described in the next section.

Results
The products of the analysis show the impassable parts of the tested terrains in multiple configurations for each vehicle. The results are computed for each direction of movement (north, east, northeast, and northwest), for two kinds of data (from LIDAR and from UAV) and for each resolution of DEM (0.5, 1, 2.5, and 5 m). Due to the large number of results, some of them are only presented briefly ( Figure 11). The arrow indicates the direction of the vehicle movement.  As far as the flat terrain is concerned (Figure 2c), due to the expectations, there was no impassable terrain indicated by the algorithm. This shows that the script works correctly.
The execution times of the script are also crucial factors of the analysis. They depend mainly on the pixel size of the model. The higher the resolution of the model is, the more time it will take to complete the execution of the script. The results of the performance analysis for tested areas (50 × 50 m) are listed in Table 4.

Discussion
To analyze the differences between the generated DEMs, in the first step, their statistics have been calculated. These are DEMs for each kind of analyzed terrain (bluff, stream, and flat) and generated on the basis of the LIDAR and UAV data that were described in Section 2. Their statistics are presented in Table 5. We can conclude from them that the differences between models from disparate sources were not high. They did not exceed 1 m value in any case. At this point, we are not able to determine which source of data is more credible. This will be verified in future research by taking field measurements. The change of the pixel size did not strongly influence the average height. There were no differences between the 0.5 and 1 m resolution models, and between the 1, 2.5, and 5 m, there was 20 cm maximally. As far as the standard deviation is concerned, its value depends on the kind of terrain, which the elevation model represents. The highest value was in the bluff terrain and the lowest was the flat terrain. In the second step of the comparison, visualizations of differences between models obtained using two sources of data were generated (the results are presented in Figure 12).
The obtained visualizations confirm the results presented in Table 5. The maximal differences between models did not exceed 0.8 m. However, the spatial distribution of disparities showed that the biggest deviations were on the areas of high height differences, which can have an impact on the results of the conducted passability analysis. Remote Sens. 2020, 12,   The obtained visualizations confirm the results presented in Table 5. The maximal differences between models did not exceed 0.8 m. However, the spatial distribution of disparities showed that the biggest deviations were on the areas of high height differences, which can have an impact on the results of the conducted passability analysis.
To make a tool that detects the microrelief objects on the digital elevation model), an accurate model must be used. There are some DEMs acquired in Poland and used widely in Polish Armed Forces and in the NATO countries: LIDAR (light detection and ranging) data, DTED (digital terrain elevation data), SRTM (shuttle radar topographic mission), and UAV (unmanned aerial vehicle) photogrammetry data. As far as the DTED and SRTM data are concerned, their horizontal accuracy is about 30 m, and the vertical accuracy is about 20 m. The pixel size of these models is about ten times larger than the size of vehicles; thus, it is worth showing that the use of these models in microrelief object detection and presented passability analysis is limited due to their relatively high resolution in comparison with the LIDAR and UAV data. This is shown in Figure 13  To make a tool that detects the microrelief objects on the digital elevation model), an accurate model must be used. There are some DEMs acquired in Poland and used widely in Polish Armed Forces and in the NATO countries: LIDAR (light detection and ranging) data, DTED (digital terrain elevation data), SRTM (shuttle radar topographic mission), and UAV (unmanned aerial vehicle) photogrammetry data. As far as the DTED and SRTM data are concerned, their horizontal accuracy is about 30 m, and the vertical accuracy is about 20 m. The pixel size of these models is about ten times larger than the size of vehicles; thus, it is worth showing that the use of these models in microrelief object detection and presented passability analysis is limited due to their relatively high resolution in comparison with the LIDAR and UAV data. This is shown in Figure 13.  The obtained visualizations confirm the results presented in Table 5. The maximal differences between models did not exceed 0.8 m. However, the spatial distribution of disparities showed that the biggest deviations were on the areas of high height differences, which can have an impact on the results of the conducted passability analysis.
To make a tool that detects the microrelief objects on the digital elevation model), an accurate model must be used. There are some DEMs acquired in Poland and used widely in Polish Armed Forces and in the NATO countries: LIDAR (light detection and ranging) data, DTED (digital terrain elevation data), SRTM (shuttle radar topographic mission), and UAV (unmanned aerial vehicle) photogrammetry data. As far as the DTED and SRTM data are concerned, their horizontal accuracy is about 30 m, and the vertical accuracy is about 20 m. The pixel size of these models is about ten times larger than the size of vehicles; thus, it is worth showing that the use of these models in microrelief object detection and presented passability analysis is limited due to their relatively high resolution in comparison with the LIDAR and UAV data. This is shown in Figure 13. The obtained results show that, in most cases, when the pixel size increased, the amount of impassable terrain indicated by the algorithm grew. There were also certain exceptions from that rule, especially when the analysis concerned the northwest and northeast direction. In this configuration, the percent of impassable terrain in the model with the 1 m pixel size was smaller than the percent in the model with the 0.5 m pixel size. This is shown on the charts in Figure 14.
The obtained results show that, in most cases, when the pixel size increased, the amount of impassable terrain indicated by the algorithm grew. There were also certain exceptions from that rule, especially when the analysis concerned the northwest and northeast direction. In this configuration, the percent of impassable terrain in the model with the 1 m pixel size was smaller than the percent in the model with the 0.5 m pixel size. This is shown on the charts in Figure 14.  From the charts (Figure 14), we can also conclude that the size of impassable terrain decreased proportionally to the increase in the ground clearance value of each vehicle. This behaviour of the results proclaims that the algorithm calculates correctly. As far as the results from the model with the 2.5 and 5 m pixel size are concerned, in most cases, there was a rapid growth of the impassable areas. Assuming that the results from the 0.5 m pixel size model are the most accurate, we can conclude that the outputs from the 2.5 and 5 m models are too generalized. This indicates that there is too much impassable terrain indicated in comparison with the results from the more accurate models. The 3D visualization of the obtained results allows us to see the described dependences (Figure 15). From the charts (Figure 14), we can also conclude that the size of impassable terrain decreased proportionally to the increase in the ground clearance value of each vehicle. This behaviour of the results proclaims that the algorithm calculates correctly. As far as the results from the model with the 2.5 and 5 m pixel size are concerned, in most cases, there was a rapid growth of the impassable areas. Assuming that the results from the 0.5 m pixel size model are the most accurate, we can conclude that the outputs from the 2.5 and 5 m models are too generalized. This indicates that there is too much impassable terrain indicated in comparison with the results from the more accurate models. The 3D visualization of the obtained results allows us to see the described dependences (Figure 15).  When the amount of the impassable terrain e.g., in two cases is the same, this does not mean that the equivalent terrain was indicated as NO GO. To check the consistency of the obtained results, an intersection analysis was conducted. The outcomes are the values of intersection between each configuration of results (from different resolution models and from different source of the data) and When the amount of the impassable terrain e.g., in two cases is the same, this does not mean that the equivalent terrain was indicated as NO GO. To check the consistency of the obtained results, an intersection analysis was conducted. The outcomes are the values of intersection between each configuration of results (from different resolution models and from different source of the data) and the values of the size of non-common impassable areas. The design of this analysis is presented in Figure 16.  When the amount of the impassable terrain e.g., in two cases is the same, this does not mean that the equivalent terrain was indicated as NO GO. To check the consistency of the obtained results, an intersection analysis was conducted. The outcomes are the values of intersection between each configuration of results (from different resolution models and from different source of the data) and the values of the size of non-common impassable areas. The design of this analysis is presented in Figure 16. Due to the large number of outputs, the discussion will concern only the results for the Honker vehicle moving northward. The values of the intersection are presented in Tables 6-11, as explained in Figure 17.  Due to the large number of outputs, the discussion will concern only the results for the Honker vehicle moving northward. The values of the intersection are presented in Tables 6-11, as explained in Figure 17.  When the amount of the impassable terrain e.g., in two cases is the same, this does not mean that the equivalent terrain was indicated as NO GO. To check the consistency of the obtained results, an intersection analysis was conducted. The outcomes are the values of intersection between each configuration of results (from different resolution models and from different source of the data) and the values of the size of non-common impassable areas. The design of this analysis is presented in Figure 16. Due to the large number of outputs, the discussion will concern only the results for the Honker vehicle moving northward. The values of the intersection are presented in Tables 6-11, as explained in Figure 17.    Table 7. Intersection between impassable areas of the Honker moving northward in bluff terrain (LIDAR data).  Table 8. Intersection between impassable areas of the Honker moving northward in bluff terrain (UAV data).   Table 10. Intersection between impassable areas of the Honker moving northward in stream terrain (LIDAR data).  Table 11. Intersection between impassable areas of the Honker moving northward in stream terrain (UAV data). From the results presented in the above tables, we can conclude that, in most cases, there was a distinct advantage of common impassable area over uncommon terrain for the compared models. However, in some cases, e.g., in Table 9, in comparison of two 0.5 m models, we can see that the intersection value is the same or smaller than the uncommon area. All differences between common and uncommon parts of the impassable terrain may result from the fact that DEMs from UAV and LIDAR are not the same as was presented in Table 2. Despite such numerical differences, the visualization of the impassable areas shows that the character of the indicated areas is similar, as is shown in Figure 18.

Honker
We present examples with different parameter configurations.
Remote Sens. 2020, 12, 4146 23 of 32 intersection value is the same or smaller than the uncommon area. All differences between common and uncommon parts of the impassable terrain may result from the fact that DEMs from UAV and LIDAR are not the same as was presented in Table 2. Despite such numerical differences, the visualization of the impassable areas shows that the character of the indicated areas is similar, as is shown in Figure 18  This conclusion is important from the point of view of the operational use of maps, because, while planning a route through a given terrain, the driver is only interested in the possibility of overcoming it. The structure and exact shape of the impassable area is not crucial information. However, detailed elevation models allow the finding of narrow, passable corridors, which eventually allow the overcoming of terrain by performing maneuvers to avoid obstacles. Examples are presented in Figure 19. Arrows indicate the possible routes of vehicles and their widths represent the real vehicle widths. This conclusion is important from the point of view of the operational use of maps, because, while planning a route through a given terrain, the driver is only interested in the possibility of overcoming it. The structure and exact shape of the impassable area is not crucial information. However, detailed elevation models allow the finding of narrow, passable corridors, which eventually allow the overcoming of terrain by performing maneuvers to avoid obstacles. Examples are presented in Figure 19. Arrows indicate the possible routes of vehicles and their widths represent the real vehicle widths. This conclusion is important from the point of view of the operational use of maps, because, while planning a route through a given terrain, the driver is only interested in the possibility of overcoming it. The structure and exact shape of the impassable area is not crucial information. However, detailed elevation models allow the finding of narrow, passable corridors, which eventually allow the overcoming of terrain by performing maneuvers to avoid obstacles. Examples are presented in Figure 19. Arrows indicate the possible routes of vehicles and their widths represent the real vehicle widths. The next question is, which output is better: the results from LIDAR or UAV data? To define the word "better" in this question, we have to find an unambiguous coefficient that will show the superiority of the obtained results. A good solution for this problem is to use the indices described in the spatial autocorrelation theory. On the basis of the spatial distribution of objects, the autocorrelation analysis measures whether and to what extent they influence each other.
The best coefficients of spatial autocorrelation are Geary's and Moran's indices [59]. Both are used for the same data context and, in most applications, give very similar results. However, they differ in the scale and the meaning of the output value. Geary's index takes only positive values from 0 to 2 (0-objects are clustered, 1-objects are distributed randomly and independently, and 2objects are dispersed) while Moran's indices are values from -1 and 1 and they are much more intuitive (1-objects are clustered, 0-objects are distributed randomly and independently, and -1objects are dispersed) [59]. For this reason, in the spatial autocorrelation analysis, we used Moran's coefficient.
In the comparison of the obtained passability areas, the results from the UAV data were perceived as more consistent and clustered than the LIDAR output. In the first example in Figure 18, the passability terrains of the following parameters: terrain, stream, pixel size, 0.5 m, and vehicle: Honker are visualized. We can clearly see that there are more separate parts from LIDAR data than from UAV. Such dependencies are visible in each pair of compared results. The confirmation of these observations are the values of Moran's coefficients. They are listed in Table 12. As anticipated, we can see that, as far as the lower resolutions are concerned (0.5 and 1 m), the values are highly positive, and there is a slight advantage of UAV over LIDAR data. In greater resolutions (2.5 and 5 m), the coefficients are the same and almost equal to the maximum value of the index. This indicates that the 0.5 and 1 m outputs were highly clustered; however, the results from The next question is, which output is better: the results from LIDAR or UAV data? To define the word "better" in this question, we have to find an unambiguous coefficient that will show the superiority of the obtained results. A good solution for this problem is to use the indices described in the spatial autocorrelation theory. On the basis of the spatial distribution of objects, the autocorrelation analysis measures whether and to what extent they influence each other.
The best coefficients of spatial autocorrelation are Geary's and Moran's indices [59]. Both are used for the same data context and, in most applications, give very similar results. However, they differ in the scale and the meaning of the output value. Geary's index takes only positive values from 0 to 2 (0-objects are clustered, 1-objects are distributed randomly and independently, and 2-objects are dispersed) while Moran's indices are values from −1 and 1 and they are much more intuitive (1-objects are clustered, 0-objects are distributed randomly and independently, and −1-objects are dispersed) [59]. For this reason, in the spatial autocorrelation analysis, we used Moran's coefficient.
In the comparison of the obtained passability areas, the results from the UAV data were perceived as more consistent and clustered than the LIDAR output. In the first example in Figure 18, the passability terrains of the following parameters: terrain, stream, pixel size, 0.5 m, and vehicle: Honker are visualized. We can clearly see that there are more separate parts from LIDAR data than from UAV. Such dependencies are visible in each pair of compared results. The confirmation of these observations are the values of Moran's coefficients. They are listed in Table 12. As anticipated, we can see that, as far as the lower resolutions are concerned (0.5 and 1 m), the values are highly positive, and there is a slight advantage of UAV over LIDAR data. In greater resolutions (2.5 and 5 m), the coefficients are the same and almost equal to the maximum value of the index. This indicates that the 0.5 and 1 m outputs were highly clustered; however, the results from the UAV data were slightly more clustered. Taking the 2.5 and 5 m values into consideration, their almost maximum Moran indexes show that the impassable areas were almost fully clustered. The reason for these very high results could be the small area of the analyzed terrain, where there was not enough space for any distribution of impassable areas.
The results showed that, in most cases, when the analysis was based on the models of the same parameters but different data sources, the percentages of impassable terrain indicated by the script were similar. The biggest variances were in the 5 m pixel size models. This was mainly caused by the inaccuracy of the model in detecting the microrelief shapes and the large area that was covered by one pixel in this model (25 m 2 -1% of the whole analyzed terrain). Considering that LIDAR data can be easily obtained from the state geodetic and cartography and the creation of the model from the UAV is much more complex and time consuming, the more optimal solution is to create DEMs for the passability analysis from the LIDAR data. The chart showing dependences between the percent of indicated impassable terrain (Honker vehicle, bluff terrain, 0.5 m pixel size) from elevation models from different data sources is presented in Figure 20.
Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 32 the UAV data were slightly more clustered. Taking the 2.5 and 5 m values into consideration, their almost maximum Moran indexes show that the impassable areas were almost fully clustered. The reason for these very high results could be the small area of the analyzed terrain, where there was not enough space for any distribution of impassable areas.
The results showed that, in most cases, when the analysis was based on the models of the same parameters but different data sources, the percentages of impassable terrain indicated by the script were similar. The biggest variances were in the 5 m pixel size models. This was mainly caused by the inaccuracy of the model in detecting the microrelief shapes and the large area that was covered by one pixel in this model (25 m 2 -1% of the whole analyzed terrain). Considering that LIDAR data can be easily obtained from the state geodetic and cartography and the creation of the model from the UAV is much more complex and time consuming, the more optimal solution is to create DEMs for the passability analysis from the LIDAR data. The chart showing dependences between the percent of indicated impassable terrain (Honker vehicle, bluff terrain, 0.5 m pixel size) from elevation models from different data sources is presented in Figure 20. In Figure 21, we can see the dependence between the time of the execution of the script and the number of analyzed points. The chart was created based on Table 4. The number of analyzed points grew due to decrease in the pixel size of the elevation model. We can see from the presented chart that the increase in the points caused an exponential growth of the running time of the script. The analyzed terrain was relatively small (50 × 50 m). If the passability analysis concerned larger terrain, the whole process would last much longer. To choose the most optimal option, we must consider two main factors: the accuracy of the results and the time of obtaining them. The conclusions are described in the next section.  In Figure 21, we can see the dependence between the time of the execution of the script and the number of analyzed points. The chart was created based on Table 4. The number of analyzed points grew due to decrease in the pixel size of the elevation model. We can see from the presented chart that the increase in the points caused an exponential growth of the running time of the script. The analyzed terrain was relatively small (50 × 50 m). If the passability analysis concerned larger terrain, the whole process would last much longer. To choose the most optimal option, we must consider two main factors: the accuracy of the results and the time of obtaining them. The conclusions are described in the next section.
Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 32 the UAV data were slightly more clustered. Taking the 2.5 and 5 m values into consideration, their almost maximum Moran indexes show that the impassable areas were almost fully clustered. The reason for these very high results could be the small area of the analyzed terrain, where there was not enough space for any distribution of impassable areas.
The results showed that, in most cases, when the analysis was based on the models of the same parameters but different data sources, the percentages of impassable terrain indicated by the script were similar. The biggest variances were in the 5 m pixel size models. This was mainly caused by the inaccuracy of the model in detecting the microrelief shapes and the large area that was covered by one pixel in this model (25 m 2 -1% of the whole analyzed terrain). Considering that LIDAR data can be easily obtained from the state geodetic and cartography and the creation of the model from the UAV is much more complex and time consuming, the more optimal solution is to create DEMs for the passability analysis from the LIDAR data. The chart showing dependences between the percent of indicated impassable terrain (Honker vehicle, bluff terrain, 0.5 m pixel size) from elevation models from different data sources is presented in Figure 20. In Figure 21, we can see the dependence between the time of the execution of the script and the number of analyzed points. The chart was created based on Table 4. The number of analyzed points grew due to decrease in the pixel size of the elevation model. We can see from the presented chart that the increase in the points caused an exponential growth of the running time of the script. The analyzed terrain was relatively small (50 × 50 m). If the passability analysis concerned larger terrain, the whole process would last much longer. To choose the most optimal option, we must consider two main factors: the accuracy of the results and the time of obtaining them. The conclusions are described in the next section.  There are also some limitations in the developed method. One of them is that the algorithm does not take the resilience of shock absorbers into consideration. This is a factor that can hardly be predicted and algorithmized due to the multiples of other factors conditioning the behaviour of shock absorbers, e.g., quality of roads (on which the vehicle moves), the driver's style of driving, or the working conditions. This method cannot be used for tanks or other vehicles that have tracks, because the wheelbase in this case is unidentifiable, and consequently we can hardly predict in which place of the ground the tracks are tangent. The presented methodology also does not take soil parameters into account.

Conclusions
The conducted research demonstrated clearly that the presented methodology with the implemented authorial method of the parallel shift of the plane was able to detect microrelief objects and, consequently, impassable areas with different resolution digital elevation models. In the analysis, we used DEMs created on the basis of LIDAR and UAV photogrammetry data in different resolutions: 0.5, 1, 2.5, and 5 m as well as the parameters of three kinds of vehicles widely applicable in the Polish Armed Forces: Honker, Star, and Humvee. The differentiation of resolutions was made to answer the question: which model is the most efficient in obtaining the best results possible in the shortest time?
The obtained results show, as expected, that the amount of the impassable area depends on the ground clearance of the vehicle. If this value grows, the number of impassable cells decreases. We can see this dependence in Figure 14. This behaviour demonstrates the correctness of the developed methodology. In different resolution outcomes, in most cases, there was a rapid growth of impassable areas between the 1 and 2.5 m pixel sizes. This proves that the output from the mentioned models is too generalized in comparison with the smaller resolution models, which are considered to be more accurate. However, the results from the 2.5 m and 5 m models should not be rejected, because they have use in the analysis conducted at the higher operational level of command. Smaller resolutions (0.5 and 1 m) show more detailed characteristics of impassable areas.
An important factor that influences the evaluation is also the intersection area, which shows the degree of consistency between the obtained outputs. In Tables 6-11, we present the numerical values of intersection. We can conclude from these that the correspondence grew with the increase in the pixel size of the used models. In smaller resolutions, the proportions between the common and uncommon impassable areas fluctuated around 50%. Despite such unsatisfactory numerical values, the visualization of outcomes showed a high consistency between them. We can see in Figure 18 that, in most corresponding cases, there was a distinguishable main character of impassable area. The difference between outputs was mainly caused by the disparity of the used DEMs.
As far as the source of the data is concerned, the results from UAV photogrammetry data visually appeared to be less dispersed, what indicates the higher quality of the data. The confirmation of this convincement is the Moran index, which provides the numerical value of the spatial autocorrelation between results. The more clustered the objects are, the higher the value of the index is. The coefficients at smaller resolutions (0.5 and 1 m) were higher for the models from the UAV data, and, at higher resolutions (2.5 and 5 m), they were the same.
This proves that DEMs created from the UAV photogrammetry data were qualitatively better. However, Figure 20 presents the dependences between the number of impassable areas in the corresponding models. This shows that, in most cases, there is more or less the same quantity of impassable area in the range of the model of the same parameters but with different data sources. The LIDAR data can be easily acquired from the State Geodetic and Cartography and UAV data must be obtained from the photogrammetry record and the aero triangulation process. This makes the UAV data more difficult and time consuming to acquire.
Taking the above into consideration, both data sources should be treated in a complementary manner. On terrain where there are available LIDAR data, which can be acquired from the governmental sources, their use provides the desired results. The unquestionable advantage of these data is the ease of obtaining them, making all works related with the preparation of passability maps possible to perform as office work.
They are available almost for the whole Poland's territory, and, consequently, detailed passability maps can be developed for a huge area. On the other hand, military or crisis management operations can be realized on areas where there is no access to such high-resolution data. Then, the acquisition of image data with use of UAVs and their development to obtain detailed DEMs can be indispensable. The time of the DEM obtaining a base of photogrammetry data can be progressively shortened due to the ongoing automation of the image processing process.
An important factor is also the time of the analysis. The dependence between the time and pixel size is presented in Figure 21. We can see the exponential growth of the time when increasing the pixel size and, consequently, the number of analyzed pixels. The most rapid rise of the time is noticeable between the 1 and 0.5 m models. All analyses in this research were conducted on relatively small areas (50 × 50 m). When it comes to larger terrains, the processing times could increase rapidly. The optimalization of the script will be studied in further research. All conducted tests allowed us to determine the most optimal solution for the passability analysis. To obtain accurate results in the optimal time, the best option was to use the 1 m model from the UAV photogrammetry data. This resolution is slightly less accurate than the 0.5 m; however, the time of the analysis is considerably higher.
Finally, the conducted analyses showed that the developed tool, which examines the terrain in terms of the passability, can be used in different aspects of the terrain analysis process. As was already mentioned, the higher resolution models (2.5 and 5 m) can be used at the operational level in planning military operations. On the other hand, 0.5 and 1 m, due to their high accuracy, may be used and implemented in unmanned ground vehicles. The knowledge of obstacles located in the terrain will increase the effectiveness of actions taken by the UGV. In this article, the issue of verification of the data was not raised and will be examined in further research.