Roughness Spectra Derived from Multi-Scale LiDAR Point Clouds of a Gravel Surface: A Comparison and Sensitivity Analysis

Abstract

The roughness spectrum (i.e., the power spectral density) is a derivative of digital terrain models (DTMs) that is used as a surface roughness descriptor in many geomorphological and physical models. Although light detection and ranging (LiDAR) has become one of the main data sources for DTM calculation, it is still unknown how roughness spectra are affected when calculated from different LiDAR point clouds, or when they are processed differently. In this paper, we used three different LiDAR point clouds of a 1 m × 10 m gravel plot to derive and analyze the roughness spectra from the interpolated DTMs. The LiDAR point clouds were acquired using terrestrial laser scanning (TLS), and laser scanning from both an unmanned aerial vehicle (ULS) and an airplane (ALS). The corresponding roughness spectra are derived first as ensemble averaged periodograms and then the spectral differences are analyzed with a dB threshold that is based on the 95% confidence intervals of the periodograms. The aim is to determine scales (spatial wavelengths) over which the analyzed spectra can be used interchangeably. The results show that one TLS scan can measure the roughness spectra for wavelengths larger than 1 cm (i.e., two times its footprint size) and up to 10 m, with spectral differences less than 0.65 dB. For the same dB threshold, the ULS and TLS spectra can be used interchangeably for wavelengths larger than about 1.2 dm (i.e., five times the ULS footprint size). However, the interpolation parameters should be optimized to make the ULS spectrum more accurate at wavelengths smaller than 1 m. The plot size was, however, too small to draw particular conclusions about ALS spectra. These results show that novel ULS data has a high potential to replace TLS for roughness spectrum calculation in many applications.

1. Introduction

Light Detection and Ranging (LiDAR) is one of the main measurement techniques used nowadays for the derivation of digital terrain models (DTM) in environmental applications. LiDAR, also referred to as laser scanning, provides highly accurate and dense 3D point measurements (point clouds) that can be acquired at different measurement scales. For example, terrestrial laser scanning (TLS) can measure micro-scale surface features (<1 m) such as pebbles or soil clods [1,2,3]. Such data typically have a coverage of a few tens of meters. A coverage of, e.g., a few kilometers is possible when a TLS system is mounted on a moving vehicle or a riverboat [4,5]. These so-called mobile mapping systems also include a global navigation satellite system (GNSS) and/or an inertial measurement unit (IMU) for georeferencing the acquired points. Laser scanning performed from airplanes (ALS) constitutes another dynamic, mobile mapping system that provides an even larger coverage, but also features a lower resolution compared with TLS. However, such data are able to depict topographic signatures at the landscape scale relevant for hydrology and geomorphology [5,6,7,8,9]. Therefore, LiDAR data are extremely useful for describing natural surfaces at several spatial scales.

One way to describe natural surfaces over different spatial scales is by using the roughness spectrum [8,10,11]. This is a DTM derivative that shows the distribution of the terrain variance over a certain range of spatial frequencies (wavelengths), where each frequency represents a sinusoidal component of the DTM in the spatial domain. More precisely, the roughness spectrum is the surface’s power spectral density; i.e., the Fourier transform of the surface’s autocorrelation function. The calculation of the roughness spectra is typically based on the Fast Fourier Transform (FFT) algorithm applied to a regularly-sampled elevation profile or a gridded DTM [8,12,13]. Hedge and Masselink [14] summarized the major steps of this methodology for time series data, but it is also valid in the case of gridded DTMs [8]. As the elevation profiles are discrete and finite, a roughness spectrum has variances only in a certain span of discrete wavelengths (a particular frequency band). The longest (DC) and the smallest (Nyquist) wavelengths of this frequency band correspond exactly to the profile length and the double of the point spacing, respectively. Thus, for a LiDAR dataset, the DC and Nyquist wavelengths are defined by the spatial coverage and sampling distance, respectively.

Roughness spectra are extensively used to characterize surface roughness in many physical and environmental models. The radar backscattering models, for example, use the roughness spectrum to better define the scattering characteristics of certain natural surfaces and to predict the amount of backscatter energy from such surfaces [15]. Furthermore, roughness spectra are used to describe the roughness length of the boundary surface and to analyze its wind velocity profile for different topographic scales [16]. Similarly, in hydraulic models, the roughness spectrum and its slope are used to characterize grain roughness and local microtopography as well as to analyze the flow velocity spectrum [17,18]. Roughness spectra are also employed to quantify the influence of the Earth’s surface processes on landscape shapes [8], or to describe the seafloor topography [19].

For most of the above applications, roughness spectra are, or can be, derived from LiDAR point clouds. These LiDAR point clouds have to be interpolated into gridded DTMs first in order to apply the Fourier-based methodology and to derive the roughness spectra. It is, however, not known which frequencies of the resulting spectra are affected by the selected DTM interpolation method. Furthermore, multi-scale LiDAR datasets have different resolutions and coverages, and thus, provide spectra in different frequency bands. To combine such complementary data, it is important to know over which spatial scales the corresponding roughness spectra can be used interchangeably. Such an analysis is important to understand better how to measure surfaces with LiDAR and how to better combine these complementary LiDAR data to get the roughness spectra over a larger frequency band.

In this study, we sampled a 1 m × 10 m gravel plot with three different LiDAR techniques and analyzed and compared the corresponding roughness spectra. The LiDAR techniques used here comprise terrestrial laser scanning (TLS), laser scanning from an unmanned aerial vehicle (ULS), and laser scanning from an airplane (ALS). An unmanned aerial vehicle (UAV) is also referred to as a remotely piloted aircraft system (RPAS), or a drone. The gravel plot was extensively sampled with TLS (2.5 h of effective scanning with 14 high-resolution scans), whereas the ALS and ULS data were acquired in the course of a topographic survey campaign with the specific aim of covering a larger river section. Therefore, the focus of the analysis was more on the TLS data. In addition to these datasets, a set of high-resolution handheld images of the gravel plot was acquired and processed in order to have a non-LiDAR-based roughness spectrum for comparison. This dataset also served as a reference for the high-frequency surface components, as the ground sampling distance (GSD) of the images was 10 times smaller than the TLS footprint size. It is noted that the ULS data were acquired with a novel RIEGL VUX-1 UAV system that was applied in a surface roughness application for the first time.

Related Work and Objectives

In this study, roughness spectra and their confidence intervals are derived by applying the Fourier-based methodology standardized by Hedge and Masselink [14]. Previous studies focused more on how particular Fourier-based steps such as profile windowing or detrending affect the shape of the roughness spectra [10,13]. Snapir et al. [12] analyzed changes of roughness spectra of an agricultural soil surface derived from multi-temporal close-range images using structure from motion. Additionally, differences between the TLS roughness spectrum with roughness spectra derived from a triangulating laser scanner and images acquired by a UAV were analyzed in References [3,20], respectively. However, a comprehensive comparison of roughness spectra derived from several measurements techniques has not yet been done.

There exist several numerical roughness studies based on synthetic surfaces. These surfaces are typically simulated from a predefined shape of the roughness spectrum (or its autocorrelation function), applying the convolution of a random (uncorrelated) Gaussian surface (or a profile) with a moving window whose coefficients are derived from the predefined roughness spectrum [21,22]. The synthetic profiles are particularly used to understand the influence of the profile length, point spacing, and other sampling parameters on different roughness coefficients, as well as the radar backscattering coefficient [23,24]. Berti et al. [25] analyzed the performance of different roughness assessment methods for automatic landslide detection on both synthetic and real surfaces. A synthetic sinusoidal surface is also used for testing the Fourier-based calculations of the autocorrelation function in a roughness analysis of seafloor roughness [26]. Landy et al. [27] developed a numerical model to simulate TLS measurements over a random, exponentially-correlated surface, which was then used to assess the influence of the occlusion errors and other TLS sampling errors. In these numerical studies, the theoretical roughness values of the synthetic surfaces are well approximated by the roughness methodology applied. However, real data and experiments come with additional complexities, which are also relevant to consider and quantify.

This study considers exclusively real, multiscale LiDAR data. As a reference, a very high-resolution DTM dataset acquired from handheld images is used. The analysis concentrates on how certain aspects of LiDAR data processing affect the shape of the resulting roughness spectrum and then compares the multiscale LiDAR spectra. More precisely, the objectives of this study are:

• To analyze the influence of basic DTM interpolation methods on TLS, ULS, and ALS roughness spectra;
• To analyze how the number of TLS scans affects the roughness spectrum;
• To compare TLS, ULS, and ALS roughness spectra.

For the DTM interpolation analysis, TLS, ULS, and ALS point clouds are interpolated using basic interpolation methods and then the difference in the shape of the corresponding spectra is analyzed. The aim is to determine the spatial wavelengths that are not affected by the interpolation method. For the sake of simplicity, complex interpolation methods (e.g., the Kriging method) are not applied in this paper as their results depend on several parameters and impose stronger restrictions on the input data.

For analyzing the number of TLS scans, different combinations of TLS scans are selected and the corresponding roughness spectra are compared with the roughness spectrum derived using all 14 TLS scans (i.e., the reference roughness spectrum in this case). The aim is to determine if a 1 m × 10 m gravel plot can be measured with a smaller number of TLS scans, while still providing an equally accurate roughness spectrum as opposed to using all the TLS scans. An additional aim is to determine which frequency band can be covered by a single TLS scan with its roughness spectrum.

A comparison of the multiscale LiDAR spectra is done for the first time here. The aim of this analysis is to understand the differences in the TLS, ULS, and ALS spectra and to determine the spatial wavelengths over which they can be used interchangeably. The spectra based on TLS and novel ULS data are separately analyzed to assess to which degree the ULS spectrum can replace the TLS spectrum.

2. Study Site and Data

Table 1. An overview of the data used in this paper. The data characteristics refer to the gravel plot area.

Technique	Sensor	Slant Range	Footprint Diameter	Mean Sampling Distance	Overlapping Footprints	Processing Stage
		(m)	(mm)	(mm)
ALS	RIEGL LSM-Q1560	640–720	~160	167	No	Georef. and block adj. strips
ULS	RIEGL VUX-1 UAV	50 2	~25	29	No	Georef. and block-adj. strips
TLS	Z+F Imager 5010c	2.5–4.5	~5	0.7	Yes	Raw data
DIM 1	Nikon D800	1.5–2.7	0.5 3	0.5 3	No	Raw data

¹ Dense Image Matching; ² Flying height above ground; ³ Ground Sampling Distance (GSD).

gives an overview of the datasets (and their basic characteristics) used in this study. All the data were collected on 26 February 2015, and refer to the same area: the 1 m × 10 m gravel plot. The ALS and ULS data were already available as georeferenced point clouds for this study. Their georeferencing is performed with a rigorous strip adjustment procedure presented in [28,29], in detail. The TLS data are processed here from raw point clouds to georeferenced scans. Therefore, the input data for the roughness spectra calculation and all further analyses include individually georeferenced ALS, ULS, and TLS point clouds.

Figure 1 shows the basic properties of the study site. The 1 m × 10 m gravel plot is located at a point bar in the tailwater region of the pre-alpine Pielach River, a tributary of the Danube in Austria (48°12′50.70″ N, 15°22′27.50″ E, WGS 84). The microtopography of the gravel plot (Figure 1d,e) mostly consists of unsorted fine (<2 cm in diameter) and medium (2–6 cm) pebbles, while several clusters of coarse (>6 cm) pebbles are also present. The plot also features microtopographic undulations at a cm to dm scale on top of a global planar slope.

2.1. ALS Data

The ALS data were acquired with a Riegl LMS-Q1560 dual-channel long-range laser scanning system which acquires two scans for each strip (x-shaped scan lines on the ground with a field of view (FoV): ±29°). The sensor additionally features a GNSS/IMU device to record the trajectory. The entire system was mounted on a Diamond DA42 aircraft. The ranging precision of this instrument is 20 mm, which is specified for a range of 250 m and estimated under the manufacturer test conditions [30].

The scanning parameters were set to cover the whole river reach and provide at least 10 points per m² for a single strip (i.e., two scans). The pulse repetition rate was 400 kHz and the flying speed was approximately 56 m/s (110 knots). The plot was covered with four strips, resulting in an average point density of 36 points per m². This corresponds to an average point sampling of approximately 17 cm. As the flying height was 640 m above the terrain, the footprint diameter was also approximately 16 cm (the beam divergence being 0.25 mrad). This meant that the footprints of the individual ALS samples were (on average) touching, or just slightly overlapping one another.

The raw ALS data had already been preprocessed (cf. [28,29]), providing a georeferenced point cloud in the global coordinate system (ETRS89/UTM33N). The ALS preprocessing involved (a) the extraction of Gaussian echoes; (b) solving the range ambiguities due to the multiple-time-around ranging of the used sensor; and (c) direct georeferencing of the extracted echoes. The preprocessing was done in the sensor manufacturer software (RiPROCESS) [29]. Additionally, a rigorous strip adjustment was done with the objective to minimize the point-to-plane distances in overlapping strip areas [28]. The standard deviation of the residuals between the individual ALS strips was reported to be 1 cm [29].

2.2. ULS Data

The ULS data were acquired with the Riegl VUX-1 UAV laser scanner mounted on the Riegl RiCOPTER, a remotely piloted octocopter. The sensor is equipped with an IMU/GNSS system to record the trajectory. The scanner itself performs online waveform processing, which provides additional attributes such as pulse shape deviation and calibrated reflectance readings for each laser echo [31]. The ranging precision of this instrument is 5 mm, which is specified for a range of 150 m and estimated under the manufacturer test conditions [32].

The scanning parameters were set to cover a larger area (~100 ha) and provide at least one point per dm² for a single ULS strip. The effective measurement rate was 350 kHz (cf. [29]) and the flying speed was approximately 8 m/s (15.5 knots). The plot was covered with ten strips (the effective FoV: 230°, cf. [29]), resulting in an average point density of 11.8 points per dm². This corresponds to an average point sampling of ~2.9 cm. As the flying height was 50 m above the terrain, the footprint diameter was ~2.5 cm (the beam divergence being 0.5 mrad). Thus, the footprints of the individual ULS samples were (on average) not overlapping with one another. As in the case of the ALS data, the raw ULS data had been preprocessed before (cf. [28,29]), providing a georeferenced point cloud in the same global coordinate system. The standard deviation of the point-to-plane residuals was reported to be below 2 cm for both the ULS strips only and the ULS and ALS strips combined [28].

2.3. TLS Data

The gravel plot was extensively scanned with a Z+F IMAGER^® 5010c mounted on a high tripod. This scanner utilizes phase shift ranging and has a small beam divergence (0.3 mrad). The beam diameter at the exit (0.1 m range) is 3.5 mm, which allows for a high-resolution scanning of close objects (maximum range < 130 m). The ranging precision of this instrument is 0.3 mm, which is specified by the manufacturer for ranges smaller than 10 m and for grey targets [33].

The scanning parameters were set to maximize the TLS sampling of the gravel plot and still have a reasonable scanning time to acquire as many scans as possible. The aim was for the footprints of the TLS samples to overlap one another in a single TLS scan. Therefore, the angular resolution was set to 0.036° (10,000 samples per full circle) and the ‘high’ scanning quality modus was selected. This TLS setting leads to a scanning time of 7 min per scan. The plot was scanned with 14 scans, 7 taken from each long side of the plot. This so-called opposite scanning minimizes occlusions due to object scan shadows. The average point density was 190.3 points per cm², which corresponds to an average point sampling of 0.7 mm. As the scanner height above the terrain was between 2.4 m and 2.6 m and the range was not larger than 4.5 m, the TLS footprint diameter was smaller than 5 mm within the plot. This meant that the TLS footprints were largely overlapping one another. For a single scan, the average sampling distance was ~2.4 mm (i.e., ~17 points per cm²) and thus, the TLS footprints are largely overlapping one another even for the single scan samples. Such a scan maximizes the resolution of the TLS data and that scanning scheme is also known as correlated scanning [34]. Furthermore, scanning from a high tripod ensured that the incidence angle was smaller than 52° for 90% of the data, which fulfilled the recommendations for soil roughness scanning [3].

The raw TLS measurements were first preprocessed and then georeferenced. The preprocessing involved the following steps: (a) the elimination of erroneous range measurements (so-called mixed pixels) caused by phase-based ranging [35]; (b) the georeferencing of the TLS scans; and (c) the export of the georeferenced points inside the gravel plot. All steps were performed using the Z+F LaserControl^® software (version 8.6) [36]. The georeferencing of the TLS scans was done indirectly, using four ground control points (GCPs) and a ‘point block adjustment’ method implemented in the above software. The GCPs were located outside the plot and scan positions, ensuring target visibility in all scans and a favorable network geometry for the co-registration. The coordinates of the GCPs were derived from the total station measurements and in the global coordinate system of the ALS and ULS data. After georeferencing the individual TLS scans, the standard deviation of the 3D distance residuals at the GCPs was 1.6 mm. The same statistic derived for individual TLS scans was 0.8 mm. The georeferenced TLS scans were then exported for an improved, global co-registration with a version of the Iterative Closest Point (ICP) algorithm [28]. The standard deviation of the point-to-plane distances between all scans after this global co-registration was 0.6 mm. This suggests a good co-registration as the ranging noise of the TLS scanner alone was ~0.3 mm [33].

2.4. Handheld Images

A set of handheld images of the plot were also acquired to calculate an additional, non-LiDAR-based, roughness spectrum. In total, 117 images were collected with the full-frame Nikon D800 camera using a 28 mm lens. The images were taken with an overlap of 70–80% from ~1.8 m height above the ground, which resulted in a GSD not larger than 0.5 mm within the plot.

The images were first oriented and then dense image matching (DIM) was applied to derive a DTM. This will be referred to as the DIM DTM. The image orientation and self-calibration were done using a bundle block adjustment method implemented in the Pix4D software [37]. The adjustment also included 18 3D GCPs uniformly distributed along the plot sides whose coordinates were measured by a total station in the global coordinate system. The resulting reprojection error was 0.13 pixels, whereas the standard deviation of the height residuals at the GCPs was 0.8 mm. Based on this camera orientation, the DIM DTM was calculated using the semi-global matching approach implemented in the SURE software [38]. We used the option for automatic interpolation of a DTM and specified its grid size to be 0.5 mm, which corresponds to the GSD of our images. This is 10 times smaller than the TLS footprint diameter, which enables the DIM DTM to serve as the reference for the high-frequency surface features.

3. Methods

3.1. Derivation of the Roughness Spectra

The derivation of the roughness spectra follows the steps for periodogram calculation recommended by Hegge and Masselink [14]. These steps involved detrending the surface heights, sampling the surface profiles (profiling), windowing (tapering), calculating one-sided periodograms, reducing the periodogram variance, and calculating confidence intervals. As our input is a set of independently georeferenced point clouds, our processing additionally involved: (a) data co-registration and (b) interpolation of point clouds into a DTM grid. Then, the DTM rows are randomly sampled to build an ensemble of gravel surface profiles to calculate the corresponding roughness spectra. The processing steps are summarized in Figure 2 and explained in the following subsections.

3.1.1. Co-Registration and Detrending

The georeferenced point clouds were given in the global coordinate system. However, a fine co-registration with a version of the Iterative Closest Point (ICP) algorithm [28] was additionally applied to eliminate any residual systematic georeferencing errors before the roughness spectra comparison. The TLS dataset served as a fixed reference for this fine co-registration and the ICP algorithm was used to optimize the transformation parameters of rigid body transformations (three rotations and three shifts) that were subsequently applied to the remaining data sets (the ULS and ALS point clouds). This processing was done using the ICP algorithm implementation in the OPALS software [39].

The point cloud detrending was done by subtracting a global planar trend from the data. Planar detrending is typically applied in practice to remove wavelengths longer than the plot length [8,12]. Detrending was also required because of the periodicity assumption in the Fourier-based analysis, i.e., to mitigate a height jump between the start-sample and end-sample of a profile [14]. The planar trend was estimated here by fitting a regression plane through the TLS points within the gravel plot. More precisely, the normal of the regression plane corresponds to the eigenvector of the smallest eigenvalue of the TLS coordinates’ covariance matrix. The detrending is thus performed by transforming each point cloud into a local object coordinate system (OCS) of the gravel plot. This (rigid-body) transformation aligned the eigenvector of the smallest eigenvalue with the z-axis of the OCS, whereas the eigenvectors of the largest and the second largest eigenvalues were aligned with the x- and y-axis of the OCS, respectively.

3.1.2. DTM Interpolation

For each point cloud in the OCS, a DTM grid was computed using the following interpolation methods: nearest neighbor (NN), triangular irregular network (TIN), and moving planes (MP). The MP interpolation has one additional parameter compared with other methods: the number of neighboring points. For example, the minimum required number of points for MP that allows for the adjustment of the plane is four, but a larger number of points can also be considered. In this study, we analyzed the sensitivity of the roughness spectra with respect to different numbers of neighbors and interpolation methods. More details about these interpolation methods can be found in Reference [40].

It should be noted that the DTM grid size was set to 0.5 mm for all DTMs in the study. This grid size is 10 times smaller than the footprint diameter of our highest resolution LiDAR data set (the TLS point cloud). Such a small grid size has certain advantages for our roughness spectra analysis. First, it enables the analysis of the behavior of the roughness spectra at frequencies around and higher than the diffraction-limit (the wavelength equal to the footprint diameter). Second, it excludes the impact of the grid size on our analysis. Furthermore, using the same grid size for all the DTMs results in the same frequency axis increments, which makes the comparison of the roughness spectra easier.

It should also be noted that our DTM heights are free from a global planar trend, as the trend was removed during the transformation into OCA (Section 3.1.1). Furthermore, the same transformation ensured that the rows and columns of the DTMs are parallel with the major and semi-major axis of the gravel plot, respectively. This means that each DTM row is 10 m long and contains 2000 samples separated from one another by 0.5 mm. These DTM rows are used subsequently to derive the roughness spectra.

3.1.3. Roughness Spectra and Periodograms

The roughness spectra were estimated as the ensemble averages of Hamming-windowed periodograms calculated from uniformly sampled DTM rows. The periodogram is the Fourier-based estimator of a roughness spectrum [8,10,12,13,41]. It is also the most commonly used one because it can be easily applied using the fast Fourier transform algorithm. Here, the periodograms are calculated using the MATLAB function periodogram [42].

Figure 3 shows the periodogram-based estimates of a roughness spectrum (derived from a TLS DTM) that are visualized as one-sided power spectral densities given in decibels (dB). Each roughness spectrum sample $n(f_{n,}{S}_n)$ is a pair consisting of a roughness spectrum value $S_n$ that refers to the spatial frequency $f_n.$ There are, in total, $N/2$+1 $(n=0,\dots ,N/2$) samples of a one-sided spectrum, where $N$ is the number of DTM row samples (and $N$ is an even number). The spectrum sample $n(f_{n,}{S}_n)$ corresponds to a surface harmonic $n$ that, in the spatial domain, is a sine wave with an amplitude $c_n$ and a spatial frequency of $f_n$. The relationship between $S_n$ and $c_n$ is given as $S_n\left[dB\right]=2·(10·{\log }_{10}{c_n}^2)/\Delta f$, where $\Delta f$ is the frequency resolution. The factor 2 is used to scale for the negative frequencies, except for the DC and Nyquist frequencies ($f_0$ and $f_{N/2}$, respectively) as they do not have a negative counterpart. The frequency resolution is given as $\Delta f=f_s/N=1/(N·\Delta x)$, where $f_s=1/\Delta x$ is the sampling frequency and $\Delta x$ is the sampling interval. In this paper, we analyze 10 m long DTM rows with a grid size of 0.5 mm, which means that $\Delta x$ = 0.5 mm and $N$ = 2000. Particular frequencies of $f_n$ are a set of the following discrete frequencies: $f_n=n/(N·\Delta x)$ for $n=0,\dots ,N/2$, where $n$ is the harmonic number. Each discrete frequency $f_n$ corresponds to the wavelength ${\lambda }_n=1/f_n$. The frequencies $f_n$ are typically plotted on a logarithmic scale (i.e., ${\log }_{10}{f}_n$) as in the bottom x-axis in Figure 3. The top x-axis shows the corresponding wavelengths. Finally, it is noted that the root-mean-square height $RMS{h}_n$ of the harmonic $n$ is $RMS{h}_n=\sqrt{2}·c_n$. In Figure 3, the right axis shows $RMS{h}_n$ values in millimetres that correspond to the roughness spectrum axis values in dB (the left axis).

3.1.4. Windowing, Variance Reduction, and Confidence Intervals

The periodogram has a large variance and spectral leakage ([43], p. 552). One way to mitigate leakage is by using the Hamming window [8,14], which is also done here. This window has a very steep roll-off at its side lobes [44] which is advantageous when analyzing the slope of a spectrum at high frequencies. The periodogram variance is reduced here in the same way as in Reference [45] by ensemble averaging. The following steps are required for computing the ensemble averaged periodogram. First, a certain number ($M$) of DTM rows is uniformly sampled to calculate the corresponding periodograms $n(f_{n,}{S}_{n,i})$, where $i=1,\dots ,M$. Then, an ensemble averaged periodogram $n(f_{n,}{\overline{S}}_n)$ is calculated by averaging individual periodogram values $S_{n,i}$ that refer to the same frequency $f_{n,}$: ${\overline{S}}_n=1/M{\displaystyle \sum }_{i=1}^M{S}_{n,i}$. There are also other approaches for variance reduction based on the binning or smoothing of the periodogram, but they reduce its resolution and are more appropriate when just a few profiles are available [14].

Figure 3 shows how ensemble averaging reduces the variance of the estimated roughness spectra. The black roughness spectrum that has the highest variance is estimated from a single-profile periodogram. On top of this periodogram, there is its 95% confidence interval (the black error bars) calculated as in [14], Equation (15):

$$S_{n,i}·\nu /{\chi }_{\nu ,1-\frac{\alpha }{2}}^2\ge S_{n,i}\ge S_{n,i}·\nu /{\chi }_{\nu ,\frac{\alpha }{2}}^2,$$

(1)

where $\alpha =0.05$ is the confidence level and $\nu$ is the degree of freedom. These confidence bounds are asymmetric, due to the chi-square distribution and the two degrees of freedom ($\upsilon =2$): $S_{n,i}~{\chi }_2^2(\upsilon ,2\upsilon )$. The red and blue roughness spectra are the ensemble averages based on the 25 and 500 periodograms, respectively. Their confidence intervals (the red and blue error bars) are more symmetrical and much smaller than that of the single periodogram case. This is because an ensemble averaged spectrum has more degrees of freedom, ${\overline{S}}_n~{\chi }_{2M}^2(2M,4M)$ compared to $S_{n,i}$ and consequently converges to the Gaussian distribution. Furthermore, this shows that the variance of the ensemble averaged spectrum is reduced compared to a single periodogram. In this paper, all the analyses are based on ensemble averaged roughness spectra using the 500 uniformly sampled periodograms.

3.2. Comparison of Roughness Spectra

To allow for a quantitative analysis, the comparison was based on the absolute differences of spectra: $\Delta S_n=abs({\overline{S}}_{n,1}-{\overline{S}}_{n,2})$. This will be referred to as the spectral difference. In a few cases, when more than two spectra are compared at once, the range statistics is used instead: $\Delta S_n=\max ({\overline{S}}_{n,j})-\min ({\overline{S}}_{n,j}),\forall j$. This will be referred to as the spectral range. Based on this $\Delta S_n$ value (also in dB as ${\overline{S}}_{n,j}$), it is then decided whether the spectra can be used interchangeably at the frequency $f_n$ (i.e., the wavelength ${\lambda }_n$). To put this decision on a more quantitative ground, $\Delta S_n$ is compared with a dB threshold that was set to the 95% confidence level of the spectra. In other words, when the $\Delta S_n$ exceeds this dB threshold, the 95% confidence intervals of the two spectra do not overlap anymore. Therefore, it is concluded that these spectra cannot be used interchangeably at the wavelength ${\lambda }_n$. It should be noted that the dB threshold is a function of only the number of profiles M (i.e., the degrees of freedom in ${\chi }_{2M}^2$) when ${\overline{S}}_n$ is analysed in the logarithmic scale (Equation (1) and Reference [14]). For our data, the dB threshold was 0.5 dB.

In total, three different comparisons are done that correspond to the three objectives set in Section 1.1. The following subsections present the spectra used in these comparisons. It should be noted that, in some of the comparisons, a spectrum derived from the DIM DTM is also used. This spectrum is referred to as the DIM spectrum.

3.2.1. Roughness Spectra and the DTM Interpolation Method

For each data set, three roughness spectra are derived using the NN-, TIN-, and four-point MP-interpolated DTMs. For the TLS data set, these three spectra are referred to as the TLS NN spectrum, the TLS TIN spectrum, and the TLS MP spectrum, respectively. The analogous naming is used for the ULS- and ALS-based spectra. For each data set, the next step is to derive the spectral range based on the corresponding NN-, TIN-, and MP-spectrum. If the spectral range exceeds the dB threshold, then the corresponding spatial wavelength (${\lambda }_{th}$) is reported and wavelengths larger than ${\lambda }_{th}$ are considered to be insensitive to the selection of the interpolation method. The NN-, TIN-, and MP-spectra are also compared with the DIM spectrum by calculating the corresponding (absolute) spectral differences. It should be noted that the influence of the interpolation method on the DIM spectrum is not analyzed because the DIM DTM was interpolated by the image matching software.

3.2.2. Roughness Spectra and the Number of TLS Scans

Our comparison analyzed both the combination of TLS scans and the individual scans. The former includes the following scan setups: a 7+7 Scan Setup, a 3+3 Scan Setup, a 2+2 Scan Setup, and a 1+1 Scan Setup. Figure 4 shows schematically the scan positions of these setups. The 7+7 Scan Setup includes a total of 14 TLS scans (seven for each of the two long sides of the plot) and its spectrum is referred to as the 7+7 Scan Setup spectrum. Spectra of the remaining setups are analogously named and they refer to six, four, and two TLS scans, respectively. The analysis of the individual TLS scans is further distinguished in two groups: the Single Centre Scans and the Single Corner Scans. The single center scans are placed approximately in the middle of a plot side (i.e., Scan 4 and Scan 14 in Figure 4). The single corner scans are placed at one of the corners of the plot (i.e., Scan 1, Scan 11, Scan 7, and Scan 17 in Figure 4).

Each of the roughness spectra corresponding to the analyzed scan setups is then compared with the 7+7 Scan Setup spectrum. The comparison was done by calculating the absolute spectral differences between these spectrum pairs. If the spectral differences are smaller than the dB threshold for all wavelengths, then the particular setup is considered to provide equally accurate roughness spectra as the 7+7 Scan Setup. In contrast, if a spectral difference exceeds the dB threshold, then the corresponding spatial wavelength (${\lambda }_{th}$) is reported.

3.2.3. Roughness Spectra and Multi-Scale LiDAR Data

The comparison of multiscale spectra is done in two steps. First, only the TLS and ULS spectra are compared as the TLS dataset has a better resolution, whereas the ULS dataset has a larger coverage. The question of how two spectra fit with one another is analyzed as in the previous cases by employing the dB threshold on the absolute spectral difference and reporting the ${\lambda }_{th}$. Then, the multi-scale spectra comparison is done for all the techniques. The LiDAR spectra (ALS, ULS, and TLS) are plotted together with the DIM spectrum and then it is analysed how they fit to one another. The observed differences and similarities of the spectrum shapes are interpreted (in the spatial domain) using color-coded DTM height maps. The LiDAR spectra used here refer to the same interpolation (the four-point MP).

For the TLS and ULS spectra, their spectral slope values (${\alpha }_{TLS}$ and ${\alpha }_{ULS}$) are additionally analysed. This parameter is the slope of the regression line fitted to a spectrum (given in the logarithmic scale) within a particular frequency band [1/${\lambda }_{n_1}$ 1/${\lambda }_{n_2}$] [13]. The spectral slopes are calculated for a series of frequency bands where ${\lambda }_{n_1}$ is fixed to the DC component (${\lambda }_{DC}$), whereas ${\lambda }_{n_2}$ is constantly decreasing up to the Nyquist wavelength ($n_2=1,\dots ,N/2$). This results in a series of TLS and ULS-based spectral slopes ${\alpha }_{TLS}({\lambda }_{n_2})$ and ${\alpha }_{ULS}({\lambda }_{n_2})$; each of them referring to a different lower wavelength ${\lambda }_{n_2}$. Finally, the relative spectral differences are calculated from the corresponding ULS and TLS spectral slopes: $\Delta \alpha ({\lambda }_{n_2})=({\alpha }_{ULS}({\lambda }_{n_2})-{\alpha }_{TLS}({\lambda }_{n_2}))/{\alpha }_{TLS}({\lambda }_{n_2})$, where $n_2=1,\dots ,N/2$. The value $\Delta \alpha ({\lambda }_{n_2})$ shows the relative spectral slope differences of TLS and ULS.

4. Results

4.1. Overview of the LiDAR DEMs and Their Quality

Before presenting the results of the roughness spectra comparison, this section gives an overview of the geometric quality of both the processed LiDAR point clouds and the interpolated DTMs.

Table 2. An overview of the geometric quality of the processed LiDAR point clouds.

	Precision				Co-Reg. Error	PointDensity
	Ranging	Angular	Planar	Vertical
	(mm)	(°)	(mm)	(mm)	(mm)	(pts/dm2)
TLS	0.3	0.007	0.4	0.5	0.6	19,030
ULS	5	0.001	2.6	5	14	11.8
ALS	20	0.001	13	20	10	0.36

presents the ranging and angular precision values that are given by the manufacturers of the used LiDAR sensors. These values are then used in the error propagation formula derived for a LiDAR polar observation that maximizes the vertical error (using the longest possible ranges from Table 1 and nadir or minimum angle observation geometry). This means that the planar and vertical precision values reported in Table 2 refer to the theoretical standard deviation values of the x, y, and z coordinates of the LiDAR points, respectively. In addition, the co-registration errors and the point density values are reported. The former refers to the standard deviation of the point-to-plane distances estimated from the flat areas.

The co-registration errors for TLS and ULS in Table 2 are larger than the theoretical vertical precision (which was expected because of error accumulation). However, the co-registration error of the ALS strips is smaller than the theoretical vertical precision of the ALS points. One reason for this is the local range averaging occurring within the much larger ALS footprint diameter compared with the ULS and the TLS footprint diameters (16 cm, 2.5 cm, and 0.5 cm, respectively).

Based on the TLS, ULS, and ALS point clouds, different interpolation methods (NN, TIN, and MP) are then applied to derive the corresponding DTMs. The relative accuracy of those LiDAR DTMs is assessed from the statistics of the vertical differences between each LiDAR DTM and the DIM DTM. These vertical differences will be referred to as the digital elevation model of differences (DoD) and

Table 3. The general statistics (mean and standard deviation) derived from the vertical residuals (DoD) between each LiDAR digital terrain model (DTM) and the DIM DTM.

	DoD NN		DoD TIN		DoD MP
	Mean	Std.	Mean	Std.	Mean	Std.
	(mm)	(mm)	(mm)	(mm)	(mm)	(mm)
TLS	–0.2	3.2	–0.2	2.8	–0.2	2.8
ULS	6.6	10.8	6.6	8.8	6.7	9.1
ALS	8.9	20.1	9.1	16.0	8.8	16.4

presents the mean and standard deviation for each interpolation method and LiDAR dataset used. The DoD values in Table 3 indicate a small bias (of a few millimeters), which is due to the data co-registration that was performed on the point clouds and not on the DTMs. However, it should be noted that the bias values are smaller than the co-registration error. It should also be noted that the observed bias does not affect further calculations as the spectral analysis is done on the detrended profiles.

All the interpolated LiDAR DTMs (nine in total: three interpolation methods × three LiDAR datasets) are presented together with the nine corresponding DoDs in the Supplementary Materials (Figure S1). A visualization of the scene is provided in Figure 5.

4.2. Sensitivity to the DTM Interpolation Method

Figure 6 compares the results when different interpolation methods (NN-, TIN-, and four-neighbor MP) are used to interpolate TLS, ULS, and ALS data. The columns from left to right in Figure 6 correspond to the three mentioned data sets, respectively. The top row shows the roughness spectra. The middle row shows spectral ranges calculated from NN-, TIN-, and MP-spectrum. The bottom row shows the absolute spectral differences of the DIM spectrum to NN-, TIN-, and MP-spectrum independently. For the simplicity of the figure, the DIM spectrum is plotted only once, in the figure with the TLS spectra (Figure 6a).

Figure 6a compares the roughness spectra that refer to the TLS data. There, the NN-, TIN-, and MP-spectrum visually agree with one another (and also with the DIM spectrum) for wavelengths larger than about 1 cm. As expected, the spectra depart from one another at smaller wavelengths. The analysis of the TLS spectral range (Figure 6b) shows that the dB threshold is actually exceeded at the wavelength of 18.3 mm. The same analyses for the ULS and ALS data (Figure 6d,e,g,h) show that their spectral ranges exceed the dB threshold much earlier, at the wavelengths of 1 m and 5 m, respectively.

Figure 6c shows the spectral differences of each TLS spectrum to the DIM spectrum. The spectral differences are the smallest for the MP TLS spectrum. The corresponding spectral differences of the ULS and ALS data are more complex, but generally, the differences of the MP spectra are the smallest. This is particularly valid around the footprint wavelength (the vertical black solid lines). It should be noted that these results refer to the MP interpolation when four neighbors are used. The analysis with a large number of neighbors was also done and the results showed that the roughness spectra departed earlier (at larger wavelengths) compared with the four-point MP. Thus, the four-point MP interpolation was used in further analysis.

4.3. Sensitivity to the Number of TLS Scans

Figure 7 shows spectra that are derived from different combinations of TLS scans as explained in Section 3.2.2. The DIM spectrum is also plotted. The results show that the setups with multiple TLS scans and the single center scans produce spectra that have a very similar shape to both the 7+7 Scan setup spectrum and the DIM spectrum (Figure 7a,b). The spectra depart (visually) from one another at wavelengths smaller than 1 cm. However, the spectral differences exceed the dB threshold at larger wavelengths. It is only the 3+3 Scan Setup spectrum for which the threshold wavelength is 3.9 mm (Figure 7a). As this wavelength is smaller than the diffraction resolution limit of the TLS data (the vertical solid line at λ = 5 mm), Figure 7 shows that the 3+3 spectrum is equally accurate as the 7+7 spectrum. Figure 7a also shows that the 2+2 and 1+1 spectra are equally accurate as the 7+7 spectra for wavelengths larger than about 4 cm (the blue and green vertical dashed lines in Figure 7a, respectively). The same is true for each of the two Single Scan Centre spectra, but only for wavelengths larger than about 5 cm. However, these four spectra exceed the threshold only locally and with the $\max (\Delta S_n)$ being 0.65 dB, which is just slightly above the dB threshold. These small differences are not visible in Figure 7a,b, and thus our reported threshold wavelengths should be considered as conservative estimates.

The shape of the Single Scan Corner spectra (Figure 7c) notably deviates from the 7+7 spectrum compared to the Single Scan Centre spectra (Figure 7b). Their threshold wavelengths are one order of magnitude larger than the ones corresponding to the other scan setups. More precisely, for wavelengths larger than 4 dm, the spectrum derived from any single scan taken from the corner of the plot is equally accurate as the 7+7 spectrum.

It should also be noted that the spectral slope derived from the corner scans already deviates for wavelengths of a few decimeters and smaller. For the central scans and other setups, the spectral slope is different only at sub centimeter wavelengths.

4.4. Comparison of TLS and ULS Roughness Spectra

Figure 8a shows the TLS and ULS spectra. Both spectra end at the wavelength that corresponds to the TLS and ULS footprint diameter (25 mm and 5 mm, respectively). The corresponding grey rectangles visualize the frequency band where the DTM interpolation influences the shape of the spectra (the ${\lambda }_{th}$ and ${\lambda }_{\mathrm{footprint}}$ values reported in

Table 4. Summary of the spectral range analysis based on the sensitivity of the basic interpolation methods (NN, TIN, and four-point MP). The dB threshold, ${\lambda }_{th}$ and $\Delta \mathrm{S}$ are introduced in Section 3.2. The values in the ${\lambda }_{\mathrm{footprint}}$ column are spatial wavelengths that correspond to the size of the laser footprint diameter. The $\Delta \mathrm{S}({\lambda }_{\mathrm{footprint}})$ column shows the spectral range values found at the ${\lambda }_{\mathrm{footprint}}$. The last two columns list the figures where the parameters are first shown and where they are later used, respectively.

Data	dB Threshold	${\mathit{\lambda }}_{\mathit{t}\mathit{h}}$	$\mathbf{\Delta }\mathbf{S}({\mathit{\lambda }}_{\mathit{t}\mathit{h}})$	${\mathit{\lambda }}_{\mathbf{footprint}}$	$\mathbf{\Delta }\mathbf{S}({\mathit{\lambda }}_{\mathbf{footprint}})$	Derived in	Used in
	(dB)	(cm)	(dB)	(cm)	(dB)
TLS	0.5	1.8	0.5	0.5	3.2	Figure 6b	Figure 8a and Figure 9a
ULS	0.5	100	0.5	2.5	8.5	Figure 6e	Figure 8a–c and Figure 9a
ALS	0.5	500	0.5	16	7.9	Figure 6h	Figure 9a

, Section 4.1). It can be seen that these rectangles (frequency bands) do not overlap, which means that the DTM interpolation does not affect the TLS spectrum in the ULS interpolation uncertainty band. Thus, the TLS spectrum can serve as the reference for the ULS spectrum on this frequency band. Figure 8b shows the absolute difference between the TLS and ULS spectra (the black line). The differences exceed the dB threshold at a wavelength of about 116 mm, almost five times that of the ULS footprint. The differences at wavelengths larger than this ${\lambda }_{th}$ are almost linear, as shown by their general trend (the red line representing a 20-sample running mean). The differences at wavelengths smaller than this ${\lambda }_{th}$ are more complex, but still below 0.8 dB.

The relative differences $\Delta \alpha$ between the ULS and TLS spectral slopes (Figure 8c) show similar behaviour as the spectral differences. Starting from a wavelength of 1 m, the slope difference constantly increases and reaches a maximum of −15% around ${\lambda }_{th}$. However, for wavelengths smaller than ${\lambda }_{th}$, the slope differences erroneously decrease and eventually change their sign.

4.5. Multi-Scale Spectra and DTMs

Figure 9a shows all the LiDAR spectra together with the DIM spectrum. As in Figure 8a, the LiDAR spectra end at the wavelengths that correspond to their footprint diameter (the ${\lambda }_{\mathrm{footprint}}$ values in Table 4) and the grey rectangles show the frequency bands affected by the DTM interpolation (the ${\lambda }_{th}$ values in Table 4). It is only for the TLS spectrum that the frequency band is shown differently. The vertical, dashed blue line marks the ${\lambda }_{th}$ as reported in Table 4 (1.8 cm, the influence of the interpolation method), whereas the longer wavelength of the grey rectangle (the left, vertical full blue line) refers to a wavelength where the TLS and DIM spectral differences exceed the dB threshold. This wavelength is just slightly smaller than 1 cm, about two times the size of the TLS footprint. Therefore, our DIM spectrum and MP TLS spectrum can be used interchangeably for wavelengths larger than 1 cm. For wavelengths smaller than 1 cm, the difference between the TLS and the DIM spectrum increases linearly.

Figure 9a also shows that the interpolation uncertainty bands of the ALS and ULS spectra overlap. This means that the ULS spectrum can serve as the reference for the ALS spectrum, but only in the non-overlapping frequency part (i.e., for wavelengths larger than 1 m). Around this wavelength, the ALS spectrum differs the most from the other spectra (Figure 9a). These modulations at the wavelengths around 1 m and larger are also visible in the spatial domain when the color-coded DTM heights are compared (Figure 9b–e). This visualization of the DTMs shows that the global microtopography pattern of the plot is similar for the DIM, TLS, and ULS DTMs, whereas, for the ALS DTM, this pattern is distorted.

The ULS spectrum differs the most from the TLS and DIM spectra at a different scale, namely at wavelengths of a few decimeters. Figure 10a–d renders the DTMs within the area marked by the black rectangles in Figure 9b–e to emphasize better the microtopographic components of the decimeter wavelengths. For the TLS and DIM DTMs, this microtopographic pattern appears very similar at this scale. For the ULS DTM, the pattern is slightly distorted. For the ALS DTM, however, the pattern is severely distorted at the scale of a few decimeters, which is expected as its footprint size is of the same order of magnitude. Finally, Figure 10e–g renders the DTMs from the area marked by the red rectangles in Figure 9b–e to emphasize the differences at the centimeter to millimeter wavelengths. Here, the ULS DTM shows a severely deformed microtopographic pattern compared to that of the TLS and DIM DTMs, which is again expected as the ULS footprint is larger than this scale. However, at this scale, the differences between the TLS and DIM DEMs also become visible. An example is a tiny branch (marked with the black arrow) that is still visible in the DIM DTM and in a nadir image of this area (Figure 10e,h, respectively), whereas this element is not visible in the TLS DTM (Figure 10f). These millimeter-scale differences between the TLS and DIM DEMs are also reflected in the corresponding spectra which start to depart from one another at wavelengths smaller than 1 cm (Figure 9a).

5. Discussion

5.1. Interpretation of the Results

This study analyses how the DTM interpolation, number of TLS scans, and multi-scale LiDAR data affect the shape of the roughness spectra. The aim of this paper was to quantify the wavelengths over which the roughness spectra can be used interchangeably.

5.1.1. DTM Interpolation

The three DTM interpolation methods (NN, TIN, and MP) affect the roughness spectrum only at the high frequencies (i.e., the spatial wavelengths close to the laser footprint diameter). The analyses from Section 4.2 (Figure 6) demonstrated that the ULS and ALS spectral ranges exceed the dB threshold much earlier (at wavelengths of 1 m and 5 m, respectively) compared with the TLS spectral range that exceeds at wavelengths of ~2 cm and smaller. This was expected as the ULS and ALS data are acquired at smaller measurement scales (coarser resolutions).

However, these results also show that the selection of the interpolation method is differently important for each LiDAR dataset. For the ULS data, the NN, TIN, or MP interpolations do not affect the roughness spectra at wavelengths larger than 1 m. Therefore, the optimization of the interpolation parameters is required to derive accurate ULS spectra at wavelengths smaller than 1 m. For the ALS data, the selection of interpolation method is particularly important as it affects almost all wavelengths. It is only the DC component and the first harmonic that are not affected by the interpolation method. One reason why these two surface components are not affected could be due to the fact that our plot length was too small for analyzing the ALS data. Therefore, a longer plot should be considered to determine accurately the threshold wavelength for the ALS data.

5.1.2. TLS Scan Setup

The analysis of the number of TLS scans is performed to further optimize the TLS setups for the measurement of the roughness spectrum. It was shown that a TLS scan taken from the plot corner may notably deform the shape of the roughness spectrum, whereas a center scan or two or more scans provide equally accurate roughness spectra as with all 14 TLS scans. It should be noted that a scan that is taken from the corner of a 1 m × 10 m plot can also be considered as a scan taken from the middle of a 1 m × 20 m plot. Based on this observation and the results shown in Figure 7b,c, the following statements can be drawn:

• a single TLS scan can measure the roughness spectra at wavelength scales of either 4 dm–20 m or 5 cm–10 m, with a maximum spectral difference less than 0.5 dB (the dB threshold value).
• a single TLS scan can measure the roughness spectra at wavelength scales of 1 cm–10 m (two orders of magnitude) but with a maximal difference less than 0.65 dB (the value of $\max (\Delta S_n)$ in this frequency band, Figure 7b and

  Table 5. Summary of the spectral range analyses based on the sensitivity of the number of TLS scans. The spectral differences ($\Delta \mathrm{S}$) are calculated relative to the 7+7 TLS spectrum (when all 14 scans are used). The dB threshold, ${\lambda }_{th}$ and $\Delta \mathrm{S}$ are introduced in Section 3.2. The values in the $\Delta \mathrm{S}({\lambda }_{th})$ column are $\Delta \mathrm{S}$ values found at the wavelength ${\lambda }_{th}$. The values in the column are ${\lambda }_{\Delta \mathrm{S}>0.65\mathrm{dB}}$ spatial wavelengths where $\Delta \mathrm{S}$ exceeds 0.65 dB. The latter dB value is the rounded $\max \left[\Delta \mathrm{S}({\lambda }_{th})\right]$, which is introduced to show that the Single Scans Centre spectra, the 2+2 spectra, and the 1+1 spectra just slightly violate the dB threshold (0.5 dB), i.e., only around ${\lambda }_{th}$ wavelengths.

  Setup	Number of Scans	dB Threshold	${\mathit{\lambda }}_{\mathit{t}\mathit{h}}$	$\Delta \mathbf{S}({\lambda }_{\mathit{t}\mathit{h}})$	${\mathit{\lambda }}_{\Delta \mathbf{S}>0.65\mathbf{dB}}$	Derived in
  		(dB)	(cm)	(dB)	(cm)
  3+3 Scan setup	6	0.5	0.4	0.5	0.3	Figure 6 and Figure 7a
  2+2 Scan Setup	4	0.5	3.6	0.57	0.5	Figure 6 and Figure 7a
  1+1 Scan Setup	2	0.5	4.1	0.62	0.8	Figure 6 and Figure 7a
  Single Scan Centre 1	1	0.5	5.3	0.61	0.6	Figure 7b
  Single Scan Centre 2	1	0.5	4.5	0.58	0.5	Figure 7b

  ).

It is noted that the 0.65 dB value of $\max (\Delta {\mathrm{S}}_n)$ is just slightly violating the threshold based on the 95% confidence interval. In addition, the above statements are valid when TLS is applied from high tripods and from the center of the plot. It is expected that a roughness spectrum derived from TLS using classical geodetic tripods is more sensitive because of the larger occlusion effects. However, this is to be analyzed in further studies.

5.1.3. Multi-Scale Data and Spectral Differences

The differences between the TLS and ULS roughness spectra that are presented in Figure 8b are linear at larger wavelengths (slightly smaller than the wavelength of 1 m) and then become more complex at smaller wavelengths (slightly larger than the ULS footprint diameter). Similar conclusions can be drawn from the relative differences $\Delta \alpha$ between the ULS and TLS spectral slopes (Figure 8c). The slope differences linearly increase and reach a maximum around ${\lambda }_{th}$. However, for wavelengths smaller than ${\lambda }_{th}$, the slope differences become more complex and eventually change their sign. On the other hand, the differences between the TLS and DIM roughness spectra, presented in Figure 8a and Figure 9, are only linear. Therefore, one reason for the different behavior of the ULS spectrum for wavelengths larger and smaller than ${\lambda }_{th}$ is that the ULS data have less overlapping footprints.

It is also noted that the ULS spectral slope values (in Figure 8c) overestimate the TLS spectral slope by about 15%, although the spectral differences are below 0.5 dB. This means that spectral slope values are more sensitive than the spectral differences. Therefore, the spectral slope may result in stronger measurement requirements compared with the requirements derived from the spectral differences.

5.2. Spectral Analysis

Investigating the roughness spectra is a method that can be used to analyze DTMs. It gives consistent results for the TLS data in the analysis of the scan setup and quantifies the loss of information when reducing the number of scans. Additionally, in the extreme case of using only one scan on the edge of the plot, the loss of quality is quantified. Likewise, the wavelength at which the spectra deviate significantly from each other due to different interpolation methods provides insight on the relevance of the choice of a certain method. The analysis of the spectra for different interpolation methods from LiDAR data at all investigated scales also allowed the conclusion that the four-point moving plane interpolation, which involves smoothing, is preferable over the NN and TIN interpolations. Therefore, the analysis of the difference to a reference spectrum provides a measure that concentrates on the resolution and damping of amplitudes of a specific DTM. This is additional information compared to other methods of DTM analysis such as the DEM of Differences (DoD), or the measures of the distribution of vertical height differences (mean, standard deviation, etc.).

5.3. The dB Threshold

In this paper, the dB threshold that corresponds to the 95% confidence interval of the ensemble averaged periodogram is used to determine the wavelengths over which the spectra can be used interchangeably. For our data, this dB threshold was 0.5 dB, which is a rather strict threshold value compared with other studies. For example, in an application of roughness spectra for microwave backscatter models, a spectral difference threshold of 2 dB was acceptable [24]. Additionally, in a seafloor roughness application, the spectral differences were analyzed using a threshold of 1 dB [19]. Therefore, the ${\lambda }_{th}$ values reported in this paper should be considered as rather conservative estimates.

A dB threshold that is tailored to a particular application would exactly define the requirements for the LiDAR measurements and the roughness spectrum calculation, but only for that particular application. A threshold based on the confidence interval of the roughness spectrum has the advantage of being application independent and furthermore, accounts for the variance of the ensemble averaged spectrum.

5.4. ALS and ULS Data

The comparison of the multi-scale LiDAR data is always challenging as such measurements involve many setup parameters. Our TLS data, for example, were optimal for the analysis because the plot was oversampled by a large number of scans and we scanned in the correlated sampling mode (i.e., with highly overlapping footprints within a single scan). This then allowed for optimizing the TLS setup for the roughness spectrum calculation. On the other hand, ALS and ULS did not offer an as extensive sampling of the plot as in the case of TLS. Therefore, to optimize the ALS and ULS setups for roughness spectra calculation, it would be important to acquire ALS and ULS data that have both highly overlapping footprints within the scan line and a large number of scans (strips) that cover the plot.

We posit the hypothesis that an extensive sampling of the plot with ALS and ULS would certainly result in more optimistic ${\lambda }_{th}$ values than the ones reported here for the ULS and ALS data. For example, in Figure 6, the spectral range value $\Delta S_{fp}$ at the footprint wavelength is almost 8 dB for the ALS and ULS data, whereas, for the TLS data, this value is about 3 dB. This shows that the ULS and ALS spectra are more sensitive to the interpolation method than the TLS spectra. One reason for this could be the lower footprint overlap of the ULS and ALS data compared with the TLS data. It would be interesting to see, e.g., if both the $\Delta S_{fp}$ and ${\lambda }_{th}$ values would be smaller than the values reported in Section 4.1 when derived from the ULS and ALS data with highly overlapping footprints. In addition, it would be interesting to see if the complex differences observed between the TLS and ULS spectra (Figure 9b) could become linear as with the differences observed between the TLS and DIM spectra (the grey rectangle with blue edges, Figure 9a). Highly overlapping footprints correspond to an extremely dense sampling of surface heights and thus, the roughness spectra will be less affected by the interpolation method. This means that an agreement between the TLS and ULS spectra can be expected to improve even to the sub-decimeter wavelengths (e.g., when the sampling distance of the ULS points is notably smaller than the ULS footprint diameter). Bearing in mind that our dB threshold is rather strict compared to other publications, it can be concluded that the ULS data has a high potential to replace the TLS data for the roughness spectrum calculation in many applications. This should be, however, shown in further experiments.

5.5. Limitations and Suggested Further Experiments

The idea behind this work was to initiate a quality analysis of the roughness spectrum calculation from LiDAR data. The outcome provided valuable conclusions on the TLS setup and the comparison of the TLS, ULS, and ALS spectra (and the range of their interchangeability). However, the experiment also revealed that both the experimental setup and the roughness spectra analysis could be improved in further studies.

One way to extend the analysis of the roughness spectra could be by propagating them further through geophysical models and analyzing how the model predictions are sensitive to changes in the roughness spectrum. This analysis could help in defining better dB thresholds for particular applications. In addition to this, further analysis could consider more complex interpolation methods, such as the Kriging method. As Section 4.2 showed, the TLS spectrum can be used as the reference in optimizing the interpolation parameters for ULS data by minimizing the differences between the ULS- and TLS-spectrum at wavelengths smaller than 1 m. Similarly, the ULS spectrum can then be used as the reference to optimize the interpolation parameters for ALS data.

Further experiments would also require improvements to the experimental setup and data acquisition. For example, the plot size should be much longer than 10 m for analyzing the ALS spectra. A plot size of 100 m, for example, would be feasible for TLS, as it can be surveyed with 10 TLS scans (according to the conclusions from Section 4.3). The ULS and ALS data should be acquired with a point spacing (both within and between the scan lines) notably smaller than their footprint diameters (the correlated scanning [34]). This means that the plot should be scanned from a large number of overlapping strips (flight lines). Such data would allow for the optimization of the ALS and ULS setups for the roughness spectrum calculation as was done here for our TLS data.

Another variable in further experiments can be the natural surface itself. Soil surface roughness is, for example, relevant for understanding the radar backscattering signal and for modelling the runoff overland flow hydraulics [12,46]. Therefore, it would be interesting to analyze whether the interpolation effects and the multi-scale data affected the roughness spectrum of a soil surface differently compared to the gravel surface.

6. Conclusions

In this paper, we used TLS, ULS, and ALS point clouds of a 1 m × 10 m gravel plot to derive and analyze roughness spectra from interpolated DTMs. The TLS, ULS, and ALS spectra were calculated as ensemble averaged periodograms. The spectral comparison was done using the dB threshold that was based on a 95% confidence interval of the ensemble averaged periodograms. The aim was to determine the scales (spatial wavelengths) over which the spectra can be used interchangeably. Furthermore, the extensive sampling of the plot with TLS allowed us to optimize the TLS measurement setup for roughness spectra calculation.

The analysis showed that one TLS scan can be used to measure 10 m long plots and to derive roughness spectra that are reliable for wavelengths larger than 5 cm (10 times the TLS footprint size). One TLS scan can also measure 20 m long plots, but then the roughness spectrum is only reliable for wavelengths larger than 4 dm. However, one TLS scan can also measure roughness spectra over 1 cm (two times the TLS footprint size) to 10 m wavelengths, but with spectral differences reaching up to 0.65 dB.

The results also showed that a TLS setup with six scans (three each for the two longer plot sides) provides equally accurate roughness spectra as any other setup that includes more than six TLS scans. Furthermore, it was shown that the TLS spectrum (based on all the scans) is not affected by the basic interpolation methods at wavelengths larger than about 2 cm (four times the TLS footprint size). Finally, the comparison of the TLS and DIM spectra showed that they agree well with one another for wavelengths larger than 1 cm (two times the TLS footprint size). These conclusions refer to a TLS scanning performed using large tripods (scan heights of about 2.5 m) and in the correlated sampling mode (overlapping footprints within a single scan).

The comparison of the ULS and TLS spectra showed that they agreed well with one another for wavelengths larger than about 1.2 dm (about five times the ULS footprint size). At these wavelengths, the ULS spectral slope overestimates the TLS spectral slope by about 15%, although the spectral difference is below 0.5 dB. To get accurate ULS spectra at wavelengths smaller than 1 m, the optimization of the interpolation parameters is required. The plot size was too small to derive more conclusions about the ALS spectrum. However, the analysis showed that only the DC component and the first harmonic (5 m wavelength) of the ALS spectrum were not affected by the interpolation methods.

The above results refer to a dB threshold that is strict compared with the dB thresholds used in other applications (e.g., 1 dB or 2 dB in seafloor roughness and microwave backscatter modelling, respectively). Thus, the above conclusions are rather conservative and show that ULS data has high potential to replace TLS data for roughness spectra calculation in many applications.