Multifractal Terrain Generation for Evaluating Autonomous Off-Road Ground Vehicles

To generate our terrain, we use a fractal technique adopted from Dawkins et al. [13]. We assign elevation values ($z$) to our DEMs at gridded $x$ and $y$ locations using the 3D W-M function [14]. This function is defined as

where

The W-M function operates as a summation of $M$ surface layers of ridge-like structures. Thus, $M$ is the number of ridges in Eqs. 1-2. The $n$ parameter is a frequency index, with $n_{max}$ being an upper cutoff frequency. We choose $n_{max}$ to be the length of our square DEMs in pixels. The $\gamma$ parameter controls the density of surface frequencies, and $D$ is the fractal dimension. The $\phi$ values are uniform random phases among the sinusoids. These random phases allow for unique DEMs to be generated using different random seeds. The parameter $L$ is a sampling length, and $G$ is an elevation scaling coefficient.

To create more anisotropic terrain, we use an integrative approach that combines three DEMs, each with different scaling and orientations of elevation features. To achieve this, we use the W-M function and vary the number of ridges ($M$), fractal dimension ($D$), and elevation scaling coefficient ($G$). Each of the three DEMs uses an $M$ value of either $16$, $32$, or $64$, determining their spatial frequency. We refer to these as the low-, mid-, and high-frequency DEMs, respectively. To smooth the DEMs, we use a low-pass Gaussian filter. The three DEMs are all scaled near the same range and with zero mean. This scaling negates the decreased DEM scaling apparent when $D$ increases, defined in Eq. 2. The three DEMs are multiplied with the pixel-wise product to combine their features into a multifractal DEM. Figure 2 shows an example of low- (2), mid- (2), and high-frequency (2) DEMs and their combined multifractal DEM (2). The multifractal DEM’s features and roughness vary more by the spatial scale and are anisotropic compared to its comprised three DEMs, as shown in Fig. 2. To meet UE’s DEM import requirements, we convert the multifractal DEM into a grayscale, $16$-bit PNG image. Note that the four DEMs in Fig. 2 correspond to the rendered terrains in Fig. 1. An example DEM grayscale image is shown in Fig. 3. The multifractal DEMs are then rendered as a terrain in UE for AGV testing, as described in Sec. 2.2.

To test if we can control terrain roughness with our method, we vary $D$ when creating the high-frequency DEM. To limit confounding variables, and because terrain roughness is found to be interdependent on both $G$ and $D$ [11], only $D$ is adjusted in the high-frequency DEM. We modify $D$ only in the high-frequency DEM because we anticipate this to have the most impact on roughness at smaller spatial scales relevant to the size of the AGV we are testing. Adjusting $D$ has been found to affect near-isotropic roughness in terrain profiles for AGV testing in prior work [12]. We reference this prior work for selecting the W-M function parameters. However, the effects of these parameters on our approach to creating more anisotropic terrain surfaces are unknown. To assess the impact of $D$ on terrain roughness, we select $D$ values of $2.3$, $2.45$, and $2.6$. These values sample $D$ between $2$ and $3$, commonly done in prior works. This sampling range preserves the fractal properties of homogeneity and scaling. [14]. We use these values to test whether we can control the resulting roughness of the terrain. The W-M function parameters defining the low-, mid-, and high-frequency DEMs in our study are found in Table 1.

To test how $D$ affects AGV traversal difficulty, we simulate AGV traversals on our generated terrain by adapting work by Young et al. [23]. Young et al. developed a simulation environment compatible with Robot Operating System (ROS) [24]. The authors use UE which includes a simulated Clearpath Husky AGV [25]. The physics engine subsystem used in UE is NVIDIA’s PhysX [26]. PhysX accurately models vehicle-terrain interaction [27]. To plan and track straight-line paths with the AGV, we adopt a short-range trajectory planner and pure-pursuit path tracker from Refs. [28, 29].

After combining the low-, mid-, and high-frequency DEMs into a multifractal DEM, we render it in UE as a 3D terrain. The DEM is first encoded as a $1009\times 1009$ pixel PNG and then is imported into UE [30]. Once imported, we scale the imported UE landscapes down by $95\%$ in the X and Y directions. This results in a resolution of $5~{}cm$ per pixel, and $50.4\times 50.4~{}m$ terrain maps. The default terrain Z (elevation) scaling in UE of $100\%$ was found to be too large to create sections of immobilizing and smooth terrain. We experimented with the AGV in terrains generated with a $D$ value of $2.45$ to choose an appropriate terrain Z scaling. We found a Z scaling of $0.75\%$ appropriate such that the terrain contained areas that were both smooth and rough. An example 3D rendered terrain in UE, with top-down and perspective views, is shown in Figs. 3, and 3, respectively.

We use Monte Carlo simulation to see if $D$ correlates with AGV traversability. We program the AGV to traverse $20$ random A-to-B straight-line-path missions at a speed of $1~{}m/s$, in $60$ different terrain maps. Since the terrains are obstacle-free, the straight-line paths represent a naive, non-terrain-aware path planning and navigation method. To ensure the missions are practical, we restrict the AGV from starting or ending paths in unreasonably rough terrain. We classify terrain as unreasonably rough by calculating the Moore neighborhood [31, 32] gradients with a $3\times 3$ pixel mask. This mask size is appropriate given that our DEM resolution is $5~{}cm$ per pixel, and the AGV’s tire radius and width are $16.5$ and $12.5~{}cm$, respectively. The maximum absolute value of the neighborhood gradients is calculated for each DEM pixel, described as

where

In Eq. 3, $Q_{grad}$ denotes the set of gradient values in the Moore neighborhood of a center pixel P. The functions $max()$ and $abs()$ represent the maximum and absolute value functions, respectively. $P_{grad}$ represents the maximum absolute valued gradient of a target (center) pixel. The Moore neighbors of pixel $P$ share a vertex or edge with it. Equation 4 shows how the gradient set calculation is done. Here, $Z_{P}$ and $Z_{Q_{i}}$ are the DEM elevation values of the center pixel $P$ and neighbor pixel $Q_{i}$. The Euclidean distance between the two pixels is calculated and represented by the $dist()$ function. This distance equals $5$ and $7.07~{}cm$ for neighbors with a shared edge and vertex, respectively.

Equation 3 is applied to all the DEM pixels, creating a gradient map. We adapt this approach for mission selection and terrain analysis. Each pixel in the gradient map is classified as having low, semi, or high roughness. Gradient values, as calculated in Eq 4, are expressed as unitless rise-over-run values. Pixels with a maximum gradient less than or equal to $50$ are classified as having low roughness. Semi-rough classified pixels have a maximum gradient greater than $50$ and less than or equal to $140$. Pixels with a maximum gradient above $140$ are classified as having high roughness. These category thresholds were determined by deploying the AGV and traversing UE terrain locations corresponding to the DEM gradient map. We found low-roughness terrain to be smooth and easily navigable. Terrain areas classified as semi-rough were moderately rough visually, but less navigable than low-roughness terrain. We found high-roughness areas very rough visually and likely to cause the AGV to become stuck or rollover. We define these categories to quantify the areas of roughness for each DEM, which will be compared with the respective $D$ values used. An example gradient map is presented in Fig. 4. Low-, semi-, and high-roughness areas are color-coded as blue, orange, and yellow, respectively.

We apply a morphological closing filter to the gradient maps to enclose rougher areas. For this, a disk structuring element is used, with effects shown in Fig. 4. We also restrict start/goal locations not to be within five meters of the map edge. This start/goal restriction is shown as a yellow border in Fig. 4. Additionally, we enforce starting locations to be in the low-roughness (blue) category. These restrictions help prevent the AGV from rolling from the deployment location before the mission starts. The restrictions also prevent the AGV from falling off the landscape in simulation. Goal locations must be in either the low or semi-rough (orange-coded) category so the A-to-B paths are more practical.

To test terrain traversability, we conduct randomized straight-line path missions with the AGV. The start and goal poses are separated by a Euclidean distance of $35~{}m$. Given our operational area with safety borders removed is $40.4\times 40.4~{}m$, the $35~{}m$ missions can be in any orientation, given the start and goal satisfy the roughness criteria. This mission length also provides an appropriate sampling of the terrain. We select a mission starting pixel randomly from the low-roughness categorized indices. We set the starting pose by randomly selecting a heading between $0$ and 360°. The candidate goal pose, if valid, lies along the heading $35~{}m$ away. Figure 4 shows $20$ example missions with start (green) and goal (red) locations. The path to traverse is shown with a line connecting each mission start and goal. Note that most of the missions in Fig. 4 include all three roughness types.

In this work, we adjust the fractal dimension ($D$) to control roughness for our terrain generation method. We generate $20$ unique DEMs for each of the three $D$ values as described in Sec. 2.1. We examine how increasing $D$ from $2.3$ to $2.45$ and $2.6$ in the high-frequency DEM affects the roughness of the $60$ output terrain maps. The composition of roughness for each DEM is analyzed. We categorize each DEM pixel as having low, semi, or high roughness, as described in Sec. 2.3. The terrain covered by each roughness category is calculated and presented as an area percentage. To avoid bias, we calculate these percentages before applying the morphological filter. We also do not include the five-meter border in the DEM analysis since it will not be traversed. The areas of low, semi, and high roughness in each terrain are grouped by the $D$ value used in the high-frequency DEM for respective terrain generation. The terrain roughness results are presented in Sec. 3.1, using bar graphs to show the correlation between $D$ and the median area of each roughness category.

We then study how terrain traversability for the AGV is correlated with $D$. The AGV traverses $20$ random, straight-line paths in each of the $60$ terrain maps. The AGV’s success rate is used as a metric to measure the traversability of the terrain maps. We define success rate as the percentage of successful missions traversed from start to goal in each map, where

In Eq. 5, $N_{success}$ is the number of successful A-to-B traversals in a terrain. The total number of A-to-B trials, equal to $20$ in our experiments, is shown as $N_{trials}$.

Due to the rough terrain, the AGV may not complete its mission due to tipping over or getting stuck while en route. To detect failures, we monitor the AGV’s orientation and location data throughout missions. If the pitch or roll of the AGV exceeds plus or minus 75°, we consider it to have tipped over. Similarly, if the AGV fails to move more than $0.2~{}m$ within $30~{}s$, we consider it to be stuck. In either case, the mission has failed and the simulation trial is ended.

We also analyze the RMS vertical acceleration and RMS pitch and roll rates only for successful missions. These dynamics serve as a measure of the AGV’s response to different terrains encountered. We analyze successful missions to understand the impact $D$ has on AGV traversability that is not captured by success rates. Successful missions are also analyzed to remove potential bias in the vehicle dynamics experienced in failed missions such as the AGV tipping over. Vertical acceleration is calculated from the AGV’s Z (elevation) time-series position data. We differentiate the position data twice to obtain vertical acceleration. Pitch and roll rates are determined by differentiating the AGV orientation time-series data. Pitch is measured relative to the AGV’s transverse axis, while roll is measured relative to its longitudinal axis. Both pitch and roll measurements are referenced from the AGV in a level orientation. Traversal data is collected at $20~{}hz$. We use RMS metrics since we are primarily interested in comparing respective metric magnitudes, and not concerned with directions. Mission traversal times are also measured for successful missions. Traversal time is the time in seconds it takes the AGV to complete a straight-line path. Since the number of successful trials varies by the terrain map used, we aggregate the successful traversal metric results by $D$ used in terrain generation. Success rate and traversal performance metrics for successful trials are presented in Sec. 3.2.

For our AGV traversal analysis, we present results in terms of the median and inner-quartile range (IQR) with boxplots. The boxplots show IQR, where

such that $Q_{3}$ and $Q_{1}$ are the third and first quartile data values, respectively. The IQR measures the spread of the center $50\%$ of a respective metric’s data points. We use median and IQR metrics since our data tends to be skewed and not normally distributed about its mean.

The roughness of the generated DEMs is analyzed by the fractal dimension used in the high-frequency DEM of our terrain generation process. Figure 5 presents the median terrain roughness for each roughness category versus the fractal dimension used. The results of our terrain generation experiments show that as $D$ increases from $2.3$ to $2.45$ and from $2.45$ to $2.6$, the median area of low-roughness terrain decreases $13.8\%$ and $7.16\%$, respectively. On the other hand, the median area of semi-rough terrain increases $11.7\%$ and $5.63\%$ as $D$ increases from $2.3$ to $2.45$, and $2.45$ to $2.6$, respectively. Additionally, the median area of high-roughness terrain increases $1.54\%$ and $3.33\%$ as $D$ increases from $2.3$ to $2.45$ and $2.45$ to $2.6$, respectively.

These results indicate that increasing $D$ in the high-frequency DEM used in our approach leads to rougher terrain on average. This is explained by the elevation assignments in our DEMs as we adjust $D$. Looking at Eq.1, as we increase $D$, the amplitude of the summed sinusoids increases exponentially, creating a rougher surface. From Eq. 2, we see that as $D$ increases, the scaling variable $A$ is scaled down exponentially, creating a flatter surface. However, we re-scale the low-, mid-, and high-frequency DEMs near the same scale before combining them. This results in the high-frequency DEM contributing and controlling the roughness realized in the multifractal DEM.

Our terrain generation method allows us to measure AGV performance over trials with varied surface characteristics. By testing AGVs in terrains generated with various combinations of fractal dimensions, we can avoid biased results and ensure the AGV can handle changes in terrain difficulty. Our terrain generation also enables us to assess AGV performance in terrains inaccessible for full-vehicle field tests. This is a crucial step in ensuring the AGV is robust enough to perform efficiently in different scenarios. In the next section, we present the results of the AGV straight-line path traversals to understand how $D$ affects traversability.

The success rate of the AGV in each terrain map is calculated based on the $20$ straight paths attempted. The success rate statistics across all terrain maps by the fractal dimension used are displayed in Fig 6. The vertical axis shows the AGV success rate in each terrain. The horizontal axis shows the corresponding $D$ value used in the high-frequency DEM for terrain generation. The red horizontal line in each boxplot indicates the median roughness composition in each set of $20$ terrains. The lower and upper box limits represent the $25\%$ ($Q_{1}$) and $75\%$ ($Q_{3}$) percentiles, respectively. Vertical dotted lines extend to the minimum and maximum terrain compositions, excluding outliers. Outliers for all boxplots are indicated as a red ”$+$” if an upper bound

or lower bound

is exceeded. Here we see that the median success rate across terrains decreases $22.5\%$ as $D$ increases from $2.3$ to $2.45$, and decreases $25\%$ as $D$ increases from $2.45$ to $2.6$.

The results of our traversal experiments show that the success rates of the AGV are affected by our control of terrain roughness. We observe that $D$ is negatively correlated with the median success rate. This is explained by the increase in the median semi- and high-roughness terrain areas and the decrease in the median low-roughness areas. This tradeoff in roughness increases the chances of the AGV encountering terrain that impedes mission completion. The results highlight the effect of greater quantities of rougher terrain on traversability. Controlling the terrain’s roughness can aid in emulating outdoor rough terrains for AGV testing. The results also highlight the need for a terrain-aware path planner to increase the likelihood of mission success.

The boxplots shown in Fig. 7 display the traversal metrics for successful missions. They depict the RMS vertical accelerations (7), RMS pitch and roll rates (7 - 7), and traversal times (7) aggregated across all $20$ terrain maps of each $D$ group. By analyzing this data, we can observe the changes in vehicle response and traversal times as $D$ increases that are not captured by success rates. Note that of the $400$ AGV trials for each $D$ value of $2.3$, $2.45$, and $2.6$ used, there were $361$, $282$, and $225$ successful trials, respectively.

As $D$ increases from $2.3$ to $2.45$ and from $2.45$ to $2.6$ we find the following. The median RMS vertical acceleration increases incrementally $0.397~{}m/s^{2}$ ($47.4\%$) and $0.281~{}m/s^{2}$ ($22.8\%$), respectively. The median RMS pitch rate increases incrementally $5.29~{}deg/s$ ($31.7\%$) and $2.13~{}deg/s$ ($9.67\%$), respectively. The median RMS roll rate increases incrementally $3.5~{}deg/s$ ($32.4\%$) and $2.02~{}deg/s$ ($14.1\%$), respectively. The median traversal time only changed when $D$ increased from $2.3$ to $2.45$, with an increase of $1~{}s$ ($2.78\%$). Note that the number of outliers in the traversal time data is $37$, $19$, and $30$ for $D$ values of $2.3$, $2.45$, and $2.6$, respectively.

The increase in the median RMS vertical acceleration and median RMS pitch and roll rates suggests that even in successful trials, the AGV experiences rougher terrain on average as $D$ increases. This data agrees with the terrain roughness results in Sec. 3.1 which shows an increase in the average areas of semi- and high-roughness terrain and a decrease in the average area of low-roughness terrain as $D$ increases. Looking at Fig. 7, it is not surprising that the median traversal times are relatively unaffected since we are analyzing only successful trials in missions defined with straight-line paths of the same Euclidean distance. However, this shows that the rougher terrain experienced by the AGV did not substantially affect median traversal times in successful missions. The traversal time outliers represent trials where the AGV’s path was impeded more significantly by rough terrain. These outliers are cases where the AGV got temporarily stuck but continued to move such that a failure was not flagged. Another cause of the outliers is rough terrain causing the AGV to deviate from the straight paths, causing longer paths to be taken.

The results of our experiments demonstrate that varying the fractal dimension in the high-frequency DEM affects the performance of the simulated AGV. Although the AGV still completes missions as $D$ increases, the average increase in its RMS vertical accelerations and pitch and roll rates show that the AGV’s kinematics are challenged more as $D$ increases. By increasing $D$, our terrain generation approach can be used to stress test a path planner or navigation method to a degraded performance or failure. This can help identify algorithm weaknesses which can then be analyzed and corrected for when operating in dangerous real-world terrain.

Note that the results discussed in this paper are specific to a modeled Clearpath Husky AGV. However, the terrain generation method outlined in this work can be utilized for vehicles of similar size and characteristics. This method can also be applied to vehicles of different sizes with the proper consideration of terrain scaling.