The Highest Gradient Model: A New Method for Analytical Assessment of the Efficiency of LiDAR-Derived Visualization Techniques for Landform Detection and Mapping

Abstract

ALS-derived raster visualization techniques have become common in recent years, opening up new possibilities for subtle landform detection in earth sciences and archaeology, but they have also introduced confusion for users. As a consequence, the choice between these visualization techniques is still mostly supported by empirical knowledge. Some attempts have been made to compare these techniques, but there is still a lack of analytical data. This work proposes a new method, based on gradient modelling and spatial statistics, to analytically assess the efficacy of these visualization techniques. A selected panel of outstanding visualization techniques was assessed first by a classic non-analytical approach, and secondly by the proposed new analytical approach. The comparison of results showed that the latter provided more detailed and objective data, not always consistent with previous empirical knowledge. These data allowed us to characterize with precision the terrain for which each visualization technique performs best. A combination of visualization techniques based on DEM manipulation (Slope and Local Relief Model) appeared to be the best choice for normal terrain morphometry, occasionally supported by illumination techniques such as Sky-View Factor or Negative Openness as a function of terrain characteristics.

1. Introduction

Airborne LiDAR (Light Detection & Ranging) technology is fast becoming one of the main methods for producing Digital Elevation Models (DEMs) as a result of its ability to efficiently collect topographic information over a large area with high precision and speed. Consequently, LiDAR is now used for a wide range of sciences and its potential has been successfully demonstrated in archaeology, geomorphology, and earth sciences, as well as other disciplines that require precise geographic information and detection of subtle landforms. LiDAR survey is especially useful when traditional surveys (aerial survey, field survey, etc.) are difficult, time-consuming, or provide low quality results [1,2].

In geomorphology, LiDAR surveys have been applied to study alluvial areas [3,4,5,6], glacial landforms [7,8], landslides dynamics in forested areas [9,10,11], or large-scale tectonic geomorphology under dense forest [2,12]. In archaeology, unknown or poorly defined sites have been discovered or better characterized through LiDAR data, particularly where vegetation cover is dense [13,14,15,16,17,18,19,20,21,22,23,24,25,26], areas are large and difficult to access [27,28], or when there are weathered remains with poor surface characteristics [29,30].

Usually, LiDAR data are filtered and classified to extract ground points, in order to produce a Digital Terrain Model (DTM) using different interpolation methods. In a second step, Derived Models (DMs) are produced from the DTM using several techniques. A broad range of DMs exist and are used as supplementary tools to improve detection and mapping of subtle landforms with high precision. These Derived Models are usually grouped into the so-called visualization techniques (VTs) which include the “classic” types (shading models, slope models, and “color cast” models), and also newer and more specific models such as Sky-view factor [19], Local Relief Model [31] or Openness [32,33]. These mostly recent (2010–) developments in VTs have resulted in a suite of DMs mainly developed for finer detection of archaeological remains and characterization, but also for geological and geomorphological use [2].

However, this diversification of VTs is confusing for the standard user. One problem is that computing methods may sometimes be complex or obscure, making the interpretation of the resulting DM far from straightforward. A more important technical constraint is that few empirical criteria exist to select the most suitable technique for each dataset and field survey. Some working strategies and guidelines have been proposed [6,34,35], but in most cases the empirical and subjective trial-and-error method for choosing the most appropriated VT remains the best option [2].

Comparisons or combinations of the different VTs, to assess their suitability, have mainly been performed according to qualitative or quantitative (but not analytical) approaches, such as visual assessment, detection counting, profile comparison, and multiple users survey [19,33,34,36,37]. These studies have provided valuable results, but have not completely eliminated sources of subjectivity and human operator bias. For example, the choice of a parameter such as the color ramp influences the representation of the different models and therefore the subjective interpretation by the user, causing confusion in assessments of the quality and suitability of each VT for a given dataset. Specialists will also often focus on certain points of their own interest, and tend to detect what they already know and unconsciously neglect other elements. Indeed, very few studies have been made to clarify the advantages and drawbacks of the different VTs. Some authors have begun to explore a complementary process, involving analytical comparison methods [35], but an in-depth analytical comparison of these techniques is still required.

This paper aims to fulfill this need for a means to evaluate LiDAR VTs by proposing a completely new analytical approach based upon spatial statistics analysis. To achieve this objective we (i) evaluated a selected panel of VTs with non-analytical approaches to provide a basis for comparison; (ii) assessed the same panel with a new fully analytical approach, the Highest Gradient Model (HGM), to test its relevance and limits; and (iii) summarized and compared the results of both approaches to extract general conclusions concerning VTs’ performance.

2. Materials and Methods

2.1. The LiDAR Survey

We tested the effectiveness and potential of applying HGM using a LiDAR dataset acquired as a part of a research project called AYPONA [38]. AYPONA is an interdisciplinary program involving geoarchaeology, paleoenvironmental analysis, environmental archaeology, and geomatics to study the plateau of Corent (Figure 1), a major archaeological site in Auvergne (French Massif Central). This site is characterized by a series of human settlements beginning in the Neolithic period and continuing until the Middle Ages [39]. The first objective of the use of LiDAR on this site was to model with a very high resolution the subtle variations of the relief, which may be indicators of the presence of archaeological remains. The second objective was to be able to map landforms to produce a precise geomorphological map, in order to identify areas of geoarchaeological or paleoenvironmental interest. Using a LiDAR survey was all the more necessary since part of the plateau and its slopes are obscured by dense forest cover.

LiDAR data were acquired by airborne laser scanner in March 2014 with a minimum of 18 emitted laser pulses per square meter. The flight altitude was 500 meters and an area of 22 km² was covered by 19 lines in two perpendicular directions. The initial LiDAR point-cloud classification was performed by the data provider to generate a raster DTM. This resulting DTM was not entirely satisfactory and we noted several classification errors such as erratic outliers, missing parts of archaeological remains and/or bare earth. In order to enhance the quality of the point cloud classification, we used the Multiscale Curvature Classification for LIDAR Data (MCC-LiDAR) algorithm, which is a command-line tool for processing discrete-return LiDAR [38,40]. The MCC algorithm applies a 3 × 3 kernel focal mean to the initial raster and returns a vector representing the mean z value. The initial raster and this vector are compared regarding a curvature tolerance Z > x + t where Z is the initial raster, t the curvature tolerance, and x the processed vector. If Z is lower than x + t, the point is classified as a ground point. The algorithm iterates over the non-classified points until it reaches a convergence threshold. The optimal bare earth classification was obtained with a curvature tolerance (t) = 0.4 and a scale parameter (s) = 0.2.

Bare earth returns were first interpolated using the natural neighbor method and converted into a regularized raster DTM, with a 0.50 m grid in the x-y plane, and with height values in meters. A 0.50 m pixel size is considered to be sufficient to preserve all the detail of bare earth point clouds and so to detect geomorphological and archaeological features [2,30,41]. The resulting raster DTM of Corent is a significant sample of a variated landscape, with a wide range of landforms and archaeological features. Therefore, it is an ideal dataset to test various relief VTs and to develop a fully analytical assessment method based on raster computing and spatial statistics.

2.2. Processing Derived Models

Within the broad range of VTs available, we selected six of the better performing ones for the visual detection of topographic anomalies. Slope (named in this paper SLOPEVIS to avoid confusion with slope as a terrain parameter) is a good compromise between extremely easy computing and reasonable results in most terrain types. Additionally, it is not dependent upon the illumination direction (unlike shading models) and interpretation is straightforward [35]. Local Relief Model or LRM [31] is an upgrade of the classical trend removal technique, considered especially useful for light relief in flat areas [30,34,35]. Both VTs are usually classified as DEM-manipulating methods.

Sky-View Factor (SVF) is an illumination technique based on the calculation of the visible sky from each position [18,19], used in urban areas but also in geomorphological mapping and archaeological remains detection. Positive and Negative Openness (OPPOS and OPNEG) are also illumination techniques based on the degree of openness of the relief at one point, used successfully in geomorphology and archaeology [32,33]. These illumination techniques have the major advantage of being independent of illumination direction, without any relief distortion. Finally I-Factor (IFACT) is a composite index from Positive and Negative Openness used in geomorphology [2]. Each VT has shown, often on an empirical basis, advantages and drawbacks depending on the terrain type, characteristics of the detection target, etc. For model processing details, see

Table 1. Selected visualization techniques, software used, and settings.

	Software	Settings	Reference
Slope (SLOPEVIS)	ArcGIS Spatial Analyst	Standard (3 × 3 sq. kernel, degrees)	-
Local Relief Model (LRM)	Open source toolbox for ArcGIS	Circular kernel (25 m)	[42]
Sky-view Factor (SVF)	Relief Visualization Toolbox 1.1	16 search directions, Circular kernel (25m), Noise removal: medium, Vertical exaggeration: 2	[18,19,43]
Positive Openness (OPPOS)	Relief Visualization Toolbox 1.1
Negative Openness (OPNEG)	Relief Visualization Toolbox 1.1
I-Factor (IFACT)	ArcGIS Raster Calculator	I = (OPPOS-OPNEG)/2	[2]

. As the kernel size was a key parameter conditioning results [34], after empirical trials of a broad range of kernel sizes, a 25 m radius circular kernel was judged optimal for all selected VTs (except for SLOPEVIS, because a 3 × 3 square kernel is the usual setting for most users) given the relatively large size of the target landforms (e.g., ancient field patterns) [31,33]. Some examples of results are given in Figure 2.

2.3. Assessing Methods for Derived Models

The selected panel of VTs was assessed with non-analytical and analytical approaches. The non-analytical approach is based on human operator detection, first with qualitative and secondly with quantitative assessment. It is easy to perform and provides valuable results, but is also subjective. We used this approach as a preliminary step to obtain a rapid assessment of VTs performance, which served as a benchmark for the results of the analytical approach developed in this work.

The analytical approach is based on spatial statistics analysis and provides fully analytical data on VTs’ performance. It was used to obtain an objective and more detailed assessment, which can be compared with the non-analytical results in order to complete the knowledge concerning the VTs, but also to highlight the advantages of this new analytical approach.

2.3.1. Non-Analytical Approach

The non-analytical approach included two methods of assessment. A qualitative assessment was based on DMs visual examination by an experienced user, who approximatively estimated (using symbols “−” to “+++”) the suitability of each VT for each landform type detection (

Table 2. Estimated performance of VTs in selected areas. Qualitative value representing the overall user impression is given by graduated symbols (− to +++).

	Rock Slumps Scars	Flood Channels	Ancient Field Patterns
SLOPEVIS	+++	−	−
LRM	+	++	+++
SVF	+	+	+
OPPOS	++	++	+++
OPNEG	+++	+++	++
IFACT	++	+++	++

). Secondly, we undertook a more exhaustive quantitative assessment, again based on visual examination but completed by manual digitalization of detected features, to count features and their measured lengths [36]. The results were presented as total cumulated length of feature detected per technique (

Table 3. Total cumulated length (in m) detected of selected landforms under different VTs. Best and worst values for each landform type are given in bold.

	Rock Slumps Scars	Flood Channels	Field Patterns
SLOPEVIS	2423	668	447
LRM	1931	2953	2700
SVF	1393	1276	2642
OPPOS	1834	1639	2770
OPNEG	1957	3261	2044
IFACT	1886	3928	2706

).

We used these two complementary methods within three test windows containing selected target landforms and archaeological remains representative of the study area diversity (Figure 3). The first window contains ancient field patterns (subtle convexities in flat relief), the second is characterized by small flood channels (alternating concavities and convexities in flat areas), and the third targets rotational landslide scars (sharp concavities in steep slopes).

For both approaches, the display parameters (color ramp, stretch, brightness, and contrast) were modified for each visual examination in order to find the best configuration for visual detection and ensure an equal-basis assessment (Figure 2 and Figure 4). This reproduces normal operator use, whereas using the same display parameters in all cases would have introduced bias as each VT has better results with different display parameters. For maximal data consistency, all the detections were performed by the same operator, with experience in the use of the different VTs, and under the same conditions.

2.3.2. Analytical Approach

Processing the Highest Gradient Model (HGM) Method

Contrast between adjacent cells is probably the most important component of visual detection, because it is the basis of shape and pattern recognition [35]. A good VT should enhance the contrast as much as possible, depending on scale and detection targets. Therefore, we propose using a contrast-based model to assess VTs. If in previous works standard deviation (StD) had been used as a proxy for contrast between cells [35], we propose the gradient parameter as an indicator of contrast. High gradient values in short distances imply high ability of a model to show clear and precise contours of subtle features. The Highest Gradient Model (HGM) determines the cells with the highest gradient for each set of DMs, and as a consequence which ones are more suitable for visual detection. A HGM for the six chosen DMs (Table 1) was built using ArcGIS software, following the steps detailed below (Figure 4).

The six initial DMs corresponding to the VTs have different units and value ranges, making them difficult to compare (Figure 4A1). In order to make them comparable, they were standardized by stretching the values to a range of 0 to 100 (Figure 4A2) using the formula I100 = [(I − Vmin)/Vmax − Vmin] × 100, where Vmin and Vmax represent the minimum and maximum values of the raster layer, I the value for a given cell, and I100 the same value stretched to the 0–100 range.

Secondly, we calculated gradients for the resulting six standardized rasters using the Spatial Analyst Slope tool with a 3 × 3 square kernel. Thirdly, we applied a trend removal (25 m circular kernel) to the six resulting gradient rasters in order to remove general trends for the benefit of local variations that were targeted by this model. From this processing we obtained six filtered gradient rasters (Figure 3A3–A5) that have the same units and similar range of values (approx. −2000 to 10,000), allowing comparisons between them.

Finally, the Highest Position tool (Spatial Analyst) was used within the methods to create the Highest Gradient Model (Figure 4B). The final HGM result appears as a six-class raster model with the same resolution as the initial DMs, and the main relief details are recognizable. Each class represents cells where one VT has higher contrast than the five others. Consequently, it is the best choice for landform detection for these cells.

Using the Highest Gradient Model for Spatial Statistics Analysis

In order to obtain detailed characterization of the terrain in which each VT performs best, we performed spatial statistics analysis combining terrain morphology layers with the HGM classes. For this we computed from the original ALS-derived DTM three raster layers representing terrain morphometry (slope, curvature, and rugosity; see Figure 4C), and overlaid them with the six classes of the HGM (one for each VT).

Slope was computed with Spatial Analyst (3 × 3 square kernel). We calculated Bolstad curvature from the original DTM with a kernel of three-cell radius circle using the “Geomorphometry and Gradient metrics” toolbox [44]. As a proxy of rugosity (a similar parameter to roughness understood as small-scale variations of amplitude in the height of a surface) we used the Surface Relief Ratio (SRR), also called Hypsometric Integral [45], which is given for a surface by:

SRR= (z(mean) − z(min))/(z(max) − z(min)).

(1)

Values usually tend to be 0.5 (a straight slope surface with no rugosity) in both directions. This was computed with the same toolbox (three-cell radius circular kernel).

Additionally we computed the standard deviation raster as a proxy for the DTM noise [35], with a circular kernel of three-cell radius (Figure 4C). The noise is considered as z variations with low amplitude and low wavelength not corresponding to the real terrain topography or judged insignificant at the scale of this work.

In order to perform a general spatial statistical analysis we first used Zonal Statistics as a Table tool (Spatial Analyst) to calculate mean and StD values of the four terrain rasters for each HGM class (Figure 4D1). When necessary, values were simplified to improve reading and interpretation (

Table 4. Highest Gradient Model zonal statistics for the panel of visualizations selected. Mean noise is computed as the StD of the StD layer. Mean curvature (Bolstad variant) is multiplied by 10³; Surface Relief Ratio value is given as 0–1 value. Values under and above the overall mean are highlighted in blue and red, respectively. Curvature and SRR values have dimensionless units.

	Noise (cm)		Slope (°)		Curvature (Bolstad) × 103		Surface Relief Ratio (SRR)
	Mean	StD	Mean	StD	Mean	StD	Mean	StD
1-SLOPEVIS	5.97	21.55	13.69	12.26	−0.27	48.39	0.484	0.12
2-LRM	1.73	12.38	7.58	7.43	−0.62	71.73	0.490	0.08
3-SVF	4.07	26.03	14.18	15.14	0.97	77.52	0.486	0.08
4-OPPOS	3.41	22.77	11.39	11.71	0.64	19.64	0.492	0.07
5-OPNEG	3.95	20.80	14.19	15.21	1.81	56.93	0.495	0.07
6-IFACT	2.65	9.52	11.69	13.58	1.58	21.05	0.497	0.07
Overall Mean	3.65	18.36	11.19	10.94	0.00	40.50	0.488	0.09

).

Secondly, for deeper and more detailed analysis we overlaid HGM classes and terrain layers using the Spatial Analyst Tabulate Area tool (Figure 4D2). For this, we first purged outliers and extreme values, and reclassified each terrain layer. Noise was reclassified in 5-mm classes, slope in 90 one-degree classes, curvature in 100 classes from −0.1 to 0.1 (−99 to +99), and SRR in 100 classes from 0 to 1 (1 to 100). Tabulate Area produced as values the total count of cells in the intersections of each terrain and HGM classes. These values were divided by the total number of cells in each terrain class in order to obtain relative values and eliminate the effect of a given terrain configuration. In the resultant graphs, the Y axis represents a relative contrast from 0 to 1 by VT for terrain classes (Figure 7, Figure 8A and Figure 9A).

Thirdly, to explore the relationship between mean noise and terrain characteristics by VT we used Zonal Statistics to overlay the values of mean noise with terrain classes, each time within a unique class of the HGM (Figure 4D3). Results were merged and shown as a distance to HGM mean noise (Figure 8B and Figure 9B). This processing was not applied to SRR due to the fact that rugosity and noise are close parameters, and therefore the analysis of noise/SRR relationship was meaningless.

3. Results

3.1. Non-Analytical Assessment Methods

3.1.1. Qualitative Image Analysis: Landform Detection

The three selected windows for visual qualitative assessing (Figure 3) included areas with anthropogenic (ancient fields patterns), fluvial (flood channels), and slope landforms (landslide scars). The results show estimated effectiveness of VTs depending on targets and terrain characteristics (Table 2). SLOPEVIS appears as a useful visualization on sharp landforms in steeper slopes as rock slump scars, but shows limited performance when used in particularly flat areas with smooth concave or convex landforms. At the opposite end of the spectrum, LRM performs extremely well with subtle landforms in flat areas, but loses effectiveness in slopes or in more marked reliefs. SVF has a homogeneous behavior in all terrains with good results. OPPOS and OPNEG work well in all terrains: the first apparently works better on convex landforms and second on concave ones. Finally, IFACT gives medium to good results in all areas, especially in landforms with alternating concavities and convexities like flood channels.

3.1.2. Counting Total Detected Landforms’ Length

Major differences appear in the total lengths of landform detected but also in the shapes, level of detail, and relative position of different segments, which could affect operator interpretations (an example applied to field patterns is given in Figure 5). The total length detected for the three selected windows is summarized in Table 3. In the case of flat areas with ancient field patterns, SLOPEVIS was once again the worst technique and OPPOS and LRM were among the best. Unexpectedly, IFACT and especially SVF showed improved behavior compared to what was previously estimated (Table 2). In areas with alternating concavities and convexities (flood channels), the best (OPNEG and IFACT) and the worst (SLOPEVIS) VTs remained unchanged. However, IFACT proved to be a little better than OPNEG, and LRM showed better results than anticipated with values close to 3000 m of feature length detected. In the steeper areas with rock slump scars, some changes in the relative performance of VTs are also noticeable: SLOPEVIS and OPNEG were the best techniques, as was previously estimated, and SVF the worst. In intermediate positions OPPOS, IFACT, and LRM appear to be much closer to OPNEG (almost equivalent in terms of total length detected) than previously estimated.

In general, the results shown in Table 3 are consistent with Table 2, but with some important differences. This shows that despite the general reliability of human assessment, the variation of display parameters for each model optimization and the human viewer can introduce important biases in technique comparisons. Additionally, differences between operators (not assessed in this work) are a well-known source of interpretation bias [35].

3.2. Analytical Assessing Methods

3.2.1. The Highest Gradient Model

The class distribution in the HGM appears to be dependent upon relief characteristics. This suggests that the best performance of a given VT is connected with precise terrain morphometry (Figure 6). The HGM deconstructs the landforms in minor homogeneous units of slope and curvature, each one corresponding with a class of the model.

The relative abundance of each class is variable. LRM has the largest number of points (very abundant in flat areas), IFACTOR very few points (largely surpassed by OPPOS and OPNEG), and SLOPEVIS, SVF, OPPOS, and OPNEG are mainly situated in the relief variations. Profiles were extracted to show the detailed functioning of the HGM and the six different VTs in different cases (Figure 6). For large and straight slopes or flat areas, LRM is almost always the best VT, whereas SLOPEVIS represents the sharp concavities and convexities with abrupt slope breaks. The short and steep slopes ranging from straight to slightly convex are represented by SVF (and sometimes also by SLOPEVIS or LRM). OPPOS and OPNEG are apparently the best VTs in cases of convexities and concavities, respectively.

3.2.2. Spatial Statistics Analysis

The spatial statistics analysis of the HGM produced quantitative results represented as statistical values and graphs. The results of Zonal Statistics analysis (mean and StD) obtained combining the six HGM classes with the three terrain morphometry and noise layers are summarized in Table 4. They give a general idea of the characteristics of terrain for which each VT obtains better results. However, a high StD, in some cases probably derived from particular distributions (e.g., not unimodal), makes interpretation difficult and suggests the need for deeper analysis. Figure 7, Figure 8A and Figure 9A provide detailed information concerning contrast/noise and contrast/terrain relationships by VT (see methods). The SRR data provided similar results to curvature; therefore, the detailed analysis of this graph is presented in Appendix A (Figure A1). The curve of total count of horizontal values was added to show the distribution of each terrain parameter. Figure 8B and Figure 9B contain complementary noise/terrain relationship data by VT (shown as deviation from mean noise). For easier interpretation, all graphs have been divided into sections as a function of curve trends and curve intersections.

General Noise Analysis

LRM and IFACT are the VTs with lowest mean noise in the HGM (1.73 and 2.65 cm) (Table 4). Both VTs also have the lowest StD (9.52 and 12.38 cm). This suggests that LRM and IFACT tend to perform better in low noise areas. SVF, OPPOS, and OPNEG appear to work better in intermediate noise areas (values around 3–4 cm) and SLOPEVIS in highest noise areas (5.97 cm). However, the high dispersion of the values (StD between 20.80 and 26.03 cm) confirms the need for detailed analysis of the noise of each VT and its contribution to contrast.

Figure 7 shows the relationship between relative contrast and noise for the six VTs, and also the noise distribution in the dataset. Most of the noise is under 50 cm (Sections I and II), and data greater than 150 cm were not considered in the analysis. Most frequent vertical noise is around 5 cm. Assuming that a VT should provide high contrast and as little noise as possible, high values in the left part of the graph (high contrast in low noise classes), and low values in the right part (low contrast in high noise classes) are considered optimal.

LRM has the best contrast for noise under 30 cm (Section I), whereas in higher noise areas (>30 cm, Sections II and III) SLOPEVIS and SVF have higher contrast values. The other VTs have low values of contrast for any noise value. This shows that LRM is the best VT concerning the contrast/noise relationship, whereas SLOPEVIS is the worst because it integrates abundant noise. SVF also shows the best contrast in Section III, suggesting that it integrates some noise. The contrast of the other three VTs is relatively independent of noise.

Terrain Parameter: Slope

LRM has the lowest mean slope (7.58°) in the HGM (Table 4). SLOPEVIS, SVF, and OPNEG have the highest mean values (13.69° to 14.19°), and OPPOS and IFACT have intermediate values between 11 and 12°. LRM has also the lowest StD (7.43°), whereas the other five VTs have moderate values (11.71° to 15.21°). This suggests that LRM is the best VT for lower slope areas, OPNEG and SVF for steeper slopes, and other VTs for intermediate slopes. The moderate to low dispersion of values also suggests well-defined contrast peaks in the given slope ranges. This was investigated through a more precise analysis of the contrast–slope relationship (Figure 8A). This graph reveals that, independently of the mean and StD values, all the VTs have different behaviors concerning slope and there is no general trend. Different VTs take the lead as the best contrast method as slope increases (successively LRM, SLOPEVIS, SVF, and OPNEG), while OPPOS and IFACT have the lowest contrasts for any slope.

In flat or low-slope areas under 15° (Section I), illumination methods (SVF, OPPOS, OPNEG, and IFACT) have low contrast values, whereas SLOPEVIS and LRM provide high contrast. LRM has higher values in these flat areas and gradually decreases with slope, whereas SLOPEVIS starts with medium contrast values and gradually increases its contrast as the slope increases. Between 15 and 30° (Section II) DEM-manipulating methods are still better, but LRM performance decreases gradually and SLOPEVIS remains more or less constant. From 30° to 60° (Section III) there is a transition area: LRM contrast falls, whereas OPNEG and especially SVF curves gradually increase with slope and SLOPEVIS remains stable. Beyond 60° (Section IV), first SVF and then OPNEG show contrast peaks, whereas the SLOPEVIS contrast decreases gradually. From 80° to 85° the curves are not interpretable due to extreme effects and outliers.

Figure 8B shows deviations from overall noise by slope class for each VT (calculated as VT specific noise–overall noise). The value 0 represents the overall mean noise for any slope class. Slopes greater than 80° were not considered in this graph due to the extreme values effect. Regardless of the slope values, SLOPEVIS always works better in noisier areas than the other techniques and with mean noise above the overall mean. This is consistent with Table 4 data, showing in addition that this trend increases with slope until 75°, but decreases in steeper slopes. By contrast, LRM is the best VT in the least noisy areas, especially for higher slopes (Section II). Under 48° (Section I) SVF, OPPOS, OPNEG, and IFACT tend to be the best VTs in areas with noise under the mean, with a trend to noise diminution. For slopes higher than 50° (Section II), these four VTs and especially OPNEG and SVF follow the opposite trend, tending to work better in areas close to the overall mean noise.

Comparing the analysis of Figure 8A,B suggests that the good performance of SLOPEVIS for medium and high slopes is always connected to higher noise, whereas the good contrast values of LRM in lower slopes corresponds to low noise values in Figure 8B. For SVF and OPNEG, good contrast values in high slopes (Figure 8A, Sections III and IV) correspond to moderate noise levels around the mean overall noise (Figure 8B Section II).

Terrain Parameter: Curvature

Mean curvature data (Table 4) show that SLOPEVIS and LRM are the only two VTs with negative mean curvatures in the HGM, suggesting that they perform better in slightly concave to straight slope areas (−0.27 and −0.62). SVF, OPPOS, OPNEG, and IFACT appear to work better in convex areas, especially OPNEG and IFACT (1.81 and 1.58). However, high StD (excepted for OPPOS and IFACT) indicate a high dispersion in the curvature values for most VTs.

This implies that mean values are not really representative and therefore precise characterization of VTs contrast as a function of curvature needs more detailed analysis (Figure 9A,B). All the VTs have different patterns concerning contrast/curvature relationship (Figure 9A). In highly concave and convex areas (Section I), SLOPEVIS contrast has much higher values than the other methods. SVF and OPNEG have much lower contrast, and for the other VTs the contrast is close to zero. For medium curvatures (Section II), there is a transition area between SLOPEVIS and LRM dominance. In areas with subtle curvatures (Section III), LRM shows the best contrast as SLOPEVIS curve continues to fall with a minimum around 0. OPNEG appears to perform slightly better in convex areas, and SVF works slightly better in concave areas. Surprisingly, OPNEG always has more contrast than OPPOS for any curvature. In general, SLOPEVIS and LRM appear to be much more significantly influenced by curvature than illumination-based techniques. Figure 9A also shows the complementarity of LRM and SLOPEVIS in terms of contrast: the first obtains better results in straight slope areas, and the latter in convex and concave areas. The SLOPEVIS curve is a good example of the benefits of this detailed analysis, because the mean value around 0 (Table 4) is clearly not representative of the true behavior of this VT concerning curvature.

Figure 9B shows deviations from overall noise by curvature class for each VT. The value 0 represents the overall mean noise for any curvature class. All VTs have a similar deviation of mean noise for slightly concave and convex areas (Section II). Independently of curvature, SLOPEVIS consistently performs better in areas with noise close to the overall mean. LRM has systematically better contrast in less noisy areas. In concave areas (Section I), OPNEG has the best contrast in areas noisier than the overall mean, whereas OPPOS and SVF have better contrast in areas with noise values under the overall mean. In convex areas (Section III) the situation is the opposite: OPNEG has better contrast in areas with noise similar to the overall mean, whereas OPPOS and SVF have the best contrast in noisier areas. For any curvatures, IFACT has intermediate noise characteristics between OPNEG and OPPOS. These data are consistent with the noise data in Table 4, providing additional details on noise variation as a function of curvature for each VT.

Comparing the analysis of Figure 9A,B shows that in straight slope areas, LRM has an optimal combination of good contrast and low noise. SLOPEVIS has the highest contrast in concave and convex areas with a noise around the overall mean noise (Figure 9B). In concave areas SVF has good balance with good contrast and medium to low noise, whereas OPNEG has less contrast and the highest noise (Figure 9B Section I). By contrast, in convex areas good performance of SVF is related to highest noise, whereas OPNEG corresponds to a noise slightly under the overall mean (Figure 9B Section III).

4. Discussion

4.1. A Detailed Assessment of ALS-Derived Raster Visualization Techniques

This work provided valuable results that allow detailed discussion pertaining to the assessment of VTs, but also methodological developments and research perspectives. As a first step, the results of the non-analytic approach highlighted differences between merely qualitative and quantitative assessing methods. Visual comparison provided only a general overview of VT performances: SLOPEVIS works better in sharp or steep areas, LRM in flat areas, and all the illumination techniques in variegated areas. Quantitative assessment globally confirmed these ideas, adding some precision mainly in relation to the performance of illumination methods. It also revealed that qualitative assessment had uncontrolled sources of bias and subjectivity in relation to display parameters and human operator survey.

The analytic approach used to assess VTs produced detailed quantitative information about each VT in relation with terrain morphology. Zonal statistical analysis (Table 4) provided a first level of analytical information. However, these data were not disaggregated enough to understand the real behavior of VTs as a function of terrain. Graphic representations of zonal tabulated data (Figure 7, Figure 8 and Figure 9 and Figure 11) allowed for a detailed characterization of the relationships between each VT and slope, curvature, rugosity, and noise. These analytical data results, summarized below, showed similarities with non-analytic assessing results but also important differences.

• SLOPEVIS has good results in medium to high slopes of any curvature except straight slopes or flat areas. As a major drawback, it integrates abundant noise independently of slope or curvature, and this noise contributes largely to the contrast in steep slopes. This excess of noise can make the detection of target features difficult or introduce bias.
• LRM shows excellent results in low slopes and flat areas, mainly in straight or low-curvature slopes. LRM integrates small amounts of noise independently of slope or curvature. The fact that the high contrast is not connected to high noise makes it an excellent technique for subtle feature detection. The main inconvenience is that its detection ability is limited to areas with predominantly low slopes or flat relief [31].
• SVF shows good results in steep slopes, in convex and especially in concave areas, where it is one of the better techniques. It also provides a reasonable performance in flat areas (sometimes better than Openness). As a major drawback, it integrates medium noise increasingly with slope and convexity, which contributes to the contrast and can make feature detection difficult.
• OPPOS has, in general, poor results compared to the other VTs, for any slope and curvature. It has only slightly better results in convex areas, where it also integrates more noise. This combination of poor contrast and noise clearly makes this VT less useful than the others.
• OPNEG showed a better than expected performance, working well in steep slopes, in concave and flat areas, but especially in convex areas, systematically outclassing OPPOS. This is consistent with the idea that these two VTs are not exactly complementary [33]. OPNEG contrast results are comparable to SVF in convex areas, but, as an additional advantage, OPNEG integrates medium to low noise in these areas, and only moderate noise in general.
• IFACT shows poor results in any terrain, probably because OPPOS and OPNEG are almost always much better. It seems to work slightly better in convex areas like OPPOS does, and to integrate medium noise. However, the very few IFACT cells in the HGM make it difficult to assess the individual behavior of this VT.

These results illustrate that each VT performs best in a terrain with precise morphometric characteristics (see Figure 10). For SLOPEVIS and LRM, results are consistent with previous investigations: this validates the HGM methodology and conforms to the analytical data ideas suggested in the literature [6,31,35].

This study has also provided new detailed data concerning illumination techniques. SVF is confirmed, as previously suggested in the literature, to work better in concavities and quite well in flat areas [19,33]. However, it also appeared to be a relatively noisy technique. OPNEG showed similar performance but with slightly better results and less noise in convex areas. OPPOS appeared as a medium to poor technique compared to others and, surprisingly, was systematically outclassed by OPNEG, even in convex areas. Finally, IFACT was the worst technique, especially when comparing with OPPOS and OPNEG results.

This work confirms that a good combination of VTs is the only solution for optimal features detection in any LiDAR dataset [33,34]. It also confirms on a fully analytical basis that this combination should be based mainly on DEM-manipulating methods such as SLOPEVIS and LRM [35,46]. This is supported by the fact that, for any dataset with a natural distribution of slopes and curvatures, they will always provide more contrast. Additionally, their noise will not be curvature-dependent, which is useful for detecting concave and convex landforms at the same time.

Considering all this, we propose precise guidelines for choosing ALS-derived raster VTs depending on terrain characteristics (Figure 10). We suggest starting with LRM and SLOPEVIS: good performance and low noise in flat areas are guaranteed by the first, and in steeper areas the latter can be successfully used. To overcome the problem of abundant noise integrated by SLOPEVIS in steeper slopes, OPNEG and/or SVF can be employed if necessary, as proposed in previous works [33]. OPNEG has good contrast and less noise in these areas and especially in convexities. As an alternative, SVF could also be used in steep slope areas implying little more noise than OPNEG, but with better results in concavities. This combination of VTs is one of the best possible choices for subtle landforms detection in most ALS datasets as it ensures an optimal contrast/noise compromise in terrain with variegated morphometry. If necessary, this general strategy for VT selection can be re-adapted to the particularities of each dataset, as a function of terrain morphometry, data quality, or target landforms.

4.2. The HGM as a Robust Method for DM Assessing: Strengths, Limitations, and Further Developments

From a methodological perspective, the HGM allowed a fully analytical quantitative assessment of ALS-derived raster VTs. The detailed developments have led to new interpretations and more precise knowledge of these techniques. Results, which are generally consistent with the existing literature and the large and variegated LiDAR dataset used in the study, validate the gradient as a proxy for contrast, and also confirm that this methodology could be applied successfully in other areas.

However, some important assumptions have to be taken into account, mainly concerning the working scale, the detection targets, and the dataset resolution. These parameters influence VTs choice but also their assessment, because most of them imply a kernel size for computing. In this work the detection targets were subtle and for large-scale landforms of any origin, implying the use of a large kernel. Therefore, the results can be extended to a large variety of common archaeological and geomorphological detection targets with similar topographic characteristics, but not systematically to small-scale landforms. On the other hand, small kernels for gradient and terrain morphology computing are adequate in most situations.

The difficulty is that the efficiency of VTs depends both on the targets’ characteristics and on the general relief of the studied area. In this study, we focused on the use of HGM to compare different VTs for a chosen kernel size. Assuming this, a controlled bias was introduced when using large kernel sizes for DM and HGM processing. It is important to note that in this study the HGM method does not resolve the problem of selecting the kernel size, but compares the results once the optimal kernel size has been found for each VT. Therefore, it is important to keep in mind that a correct application requires a conscious kernel choice for the considered VTs. However, HGM could also be used to compare, for a given VT, the influence of the kernel size on the detection efficiency, and maybe show a correlation between the kernel size and the size of the target.

The HGM method also has some important limitations to consider. The most important, the noise signal problem, is a common issue in signal treatment and analysis. In this kind of dataset, even if DMs aim to filter a large part of the noise, it is difficult to separate signal (real terrain rugosity) from small-scale noise. This implies that we will have to make some difficult decisions in filtering and data processing, in complex multi-scaled terrain morphometry parametrization, and show caution in interpreting results. In addition to this, rugosity and roughness are similar parameters but do not translate exactly into equivalent topographic realities. In any case, these two similar values are proven to be (at the scale considered in this work) of little influence regarding terrain curvature. Other questions that have not been deeply analyzed in this paper, but are of importance in the visualization of spatial data, are the settings of display parameters and the question of position mismatch and distortions introduced by different DMs. All these general issues are beyond the scope of this paper.

This HGM method also suggests the possibility of using landform signatures, understood as the complex patterns that precisely characterize each landform type in the HGM. As shown in this work, a precise pattern of slopes, concavities, and convexities has a precise signature in the HGM, which allows rapid identification of the type of landform, but also its state of weathering (e.g., well-preserved terraces shown in Figure 6B have different signatures than eroded terraces in Figure 6D). Deeper spatial statistical analysis, such as multivariate statistics, could also provide detailed information about VTs and complex microtopography interactions. This could lead to a better characterization of the “microtopographic niche” of each technique, understood as the subtle micromorphometric conditions where a VT is always best among others. These developments could have future implications for VT analysis but also for the automated detection of landforms or archaeological remains and mapping, a main challenge in LiDAR data processing [41,47].

5. Conclusions

The main objective of this work was to propose a new approach with the capacity to produce fully-analytic objective results in ALS-derived raster visualization techniques assessment. The Highest Gradient Model method was applied to a panel of six outstanding relief visualization techniques for producing reliable analytical knowledge on these techniques by comparing their functioning, within a large ALS dataset with variegated relief features. Traditional (non-analytical) qualitative and quantitative analysis of the same visualization techniques panel were also performed for methodological comparison. The spatial statistics analysis of the Highest Gradient Model provided much more detailed and objective data than the traditional approach.

A part of the previous knowledge concerning these six visualization techniques’ relative performance could be validated and consolidated. Detailed analysis also produced new information concerning recent illumination techniques such as Sky-View Factor and Openness. DEM-manipulating visualization techniques such as Slope and Local Relief Model appear to be the most powerful techniques in most cases. In specific situations, Sky-View Factor and Negative Openness can also be useful. Positive Openness and I-Factor revealed poor performance compared to the other techniques. These new results, only partly consistent with traditional non-analytical assessment results (from this work and from the literature), confirmed the importance of these kinds of analytical approaches.

As a major contribution, this work provided a robust characterization of different visualization techniques’ performance depending on the terrain characteristics, allowing us to propose objective guidelines for visualization technique choice and workflow planning when processing ALS datasets for subtle landform detection in earth sciences or archaeology.

The Highest Gradient Model method and its statistical analysis also proposed a new methodological path with future potential even if some limits remain, mainly concerning signal processing or assumptions in model computing. Advanced spatial statistical analysis of ALS-derived datasets could become a valuable tool to aid decisions in visualization and mapping processes. The microtopographic analysis and the landform signatures could also provide additional strength to automated extraction and analysis of relief features, which is probably the most important forthcoming step in the development of LiDAR data visualization.