MODELLING PERMAFROST TERRAIN USING KINEMATIC, DUAL-WAVELENGTH LASER SCANNING

In this paper we introduce the first dual-wavelength, kinematic backpack laser scanning system and its application on high resolution 3D terrain modelling of permafrost landforms. We discuss the data processing pipeline from acquisition to preparation, system calibration and terrain model process. Topographic information is vital for planning and monitoring tasks in urban planning, road construction for mass calculations, and mitigation of flood and wind related risks by structural design in coastal areas. 3D data gives possibility to understand natural processes inducing changes in the terrain, such as the cycles of thaw-freeze in permafrost regions. Through an application case on permafrost landforms in the Arctic we present the field practices and data processing applied, characterize the data output and discuss the precision and accuracy of the base station, trajectory and point cloud data. Two pulsed time of flight ranging, high performance mobile laser scanners were used in combination with a near navigation grade GNSS-IMU positioning on a kinematic backpack platform. The study shows that with a high-end system 15 mm absolute accuracy of 3D data could be achieved using PPP processing for the GNSS base station and multi-pass differential trajectory post-processing. The PPP solution shows millimetre level agreement (Easting 6 mm, Northing 4 mm, and elevation 8 mm standard deviations) for the base station coordinates over an 11 day period. The point cloud residual standard deviation for angular boresight misalignment was 27 mm. The absolute distance between ground surfaces from interactive analysis was 17 mm with 13 mm standard deviation (n=64). The proposed backpack laser scanning provides accurate and precise 3D data and performance over considerable land surface area for detailed elevation modelling and analysis of the morphology of features of interest. The high density point cloud data permits fusion of the dual-wavelength lidar reflectance data into spectral products.  Corresponding author


INTRODUCTION
Permafrost underlies 24% of the Earth's land area and is a major control for the generation of patterned ground and other terrain anomalies that can also be found on Mars (Andres et al., 2019). Periglacial features with distinctive topographic signatures, such as polygonal patterned ground and brainterrain, cover wide areas across the mid-latitudes of Mars (Balme et al., 2009;Levy et al., 2009ab;Williams et al., 2009;Zanetti et al., 2010). The origin of polygonal patterned ground and its relationship to ice on both Earth and on Mars remains debated, although a wide range of possible terrestrial analogues have been documented (Levy et al., 2009b). Patterned ground are periglacial features of interest that can provide insight into past climate, water availability, and geologic substrate on both the Earth and Mars (Washburn, 1980;Lee et al., 1998). Polygons are a common periglacial landform that are thought to be produced from the thermal contraction and subsequent cracking of frozen ground (Bramson et al., 2015). The factors that drive polygon morphology are not well understood and it is unclear what influence the subsurface ice content and geometry have on the surface expression (Hibbard et al., 2020a;French, H.M., 2017). Polygons influence the organization of water on the landscape because sloping and connected polygon troughs are flow paths for surface water (Jellinek et al., 2019).
In combination with GPR (Ground Penetrating Radar), the high resolution terrain model provides possibility to reveal the structure and faults beneath the surfacelayered soil, rocks, and ground ice. The laser scanning data allows one to correct the GPR data for topographic effects, and improve GPR data interpretation. Accurate terrain modelling allows the analysis of topographic properties and factors of each site possibly affecting the occurrence, interaction and development of such landforms induced by the micro topography in terms of slope, orientation, surface composition, freeze-thaw activity, active layer depth, and presence of melt water. Precise topography also opens possibilities for remote sensing data interpretation, and vital ground truth data for e.g. satellite radar signal processing (Zanetti et al., 2017). As we search for near-surface ice on Mars for future in-situ resource utilization, it is pertinent to understand the relationship that surface morphology has with subsurface ice content on polygon ground (Hibbard et al., 2020a) and gullies (Osinski et al., 2020).
Technological development of backpack lidar for topographic applications dates 10 years back (see e.g., Kukko et al., 2012) with early works on terrain modelling of the geomorphology of sub-arctic rivers and permafrost palsa landforms (see also Alho et al., 2009;Alho et al., 2011), and grain scale topographic mapping (Wang et al., 2013;Kukko et al., 2015). The technology and instrumentation has been continually developing during the years bringing significant reduction in sensor size and increase in the capacity and performance in terms of accuracy and measurement rates.
Topographic elevation modelling has been the primary task for mapping for long time. Such data have largely been produced by stereophotogrammetric procedures, but more recently using laser scanning, or lidar, from the air to derive the ground elevation (Pfeiffer et al., 1998;Axelsson, 2000;Vosselman, 2000;Hyyppä et al., 2000). Laser scanning is an active measurement/surveying technology to use short, nanosecond level laser pulses to probe the surface beneath (Wehr and Lohr, 1999). This requires the use of advanced GNSS-IMU technology (Global Navigation Satellite system-Inertial Measurement Unit) to observe the platform movements and reconstruct the system trajectory, i.e. position and orientation, accurately as a function of time. This allows processing of raw laser scanning data, i.e. range, scan angle (beam direction) and target reflectance, into point clouds with precise 3D locations.
Airborne laser scanning systems produce typically accuracy of 10-100 cm at relatively low point densities, typically 1-50 points per square meters. Coarse DEMs are however insufficient for analysing rugged morphology and for abundant presence of loose sediments of different sizes ranging from small scree to large boulders, hampering detailed analysis and modelling of shallow recent land surface processes, such as runoff, bioturbation or solifluction, which operate on a finer spatial scale than the coarse DEMs represent (Šašak et al., 2019). Further, the airborne applications, be it a full size ALS or a drone operated lidar, oftentimes lack perspective to reach rugged, concave shapes and other features of interest of the terrain, usually result of water flow or wind, or cracking of the ground ice, rock and soil. With ground perspective using either TLS or MLS (Terrestrial or Mobile Laser Scanning) allows such measurements, though that poses a data processing challenge for the community to address for permitting automated point classification schemes in the presence of vegetation, for example.
A recent study discusses combination of structure from motion using UAS (Unmanned Aerial System) imagery and TLS in modelling Alpine regions to overcome perspective induced occlusion issues in the data (Šašak et al., 2019). Though precise, TLS poses a major limiting challenge in covering large spatial extents: the data needs significant overlap for mutual registration, usually external targets are needed to guarantee sufficient match, and the time needed to make scans of sufficient resolution may take a significant time toll. UAS takes experienced operator, requires suitable landing sites, and complexity of UAS brings challenges to field operations, especially in the harsh Arctic climate conditions. Monitoring dynamic geomorphologic phenomena requires recurrent surveys depending on the occurrence of the event (Šašak et al., 2019). Repeatability for monitoring of e.g. erosion processes is hampered to the lack of temporal accessibility, but also lack of precision in comparison to the magnitude of the phenomena under scrutiny, and thus resolution power of subtle changes at the level of millimetres, relevant magnitude of some slow processes occurring e.g. in permafrost environments. Precision of the data from aerial perspective has been lately improved by emerging UAS sensors. The current level of precision of sensors available and suitable for limited payload UASs does not reach that of backpack and vehicle applications.
Backpack and all-terrain vehicle mounted kinematic laser scanning systems have been used before to produce e.g. forest information data for precision forestry inventory purposes (Liang et al., 2018;Wang et al., 2013). These data were collected under the forest canopy and were affected by GNSS signal deterioration (Kaartinen et al., 2012). In the Arctic polar deserts canopy occlusion is non-existent, and numerous GNSS satellites are available at all times. We investigate use of PPP processing for solving the base station coordinates.
The study area of this paper on Axel Heiberg Island, located in the Canadian High Arctic (79°22'21''N, 87°48'49''W), situates in the continuous permafrost zone where periglacial features are widespread (French, 2017). This paper presents field work practices and results of topographic modelling using Akhka-R4DW kinematic dual-wavelength laser scanning system of permafrost landforms seen abundantly in the polar deserts of the Earth. The study of these landforms provides insights into our climatic past, but may possess answers to questions regarding the climate and history of extra-terrestrial worlds in our solar system, and tools to monitor and understand the implications on the coming change in our own.

AKHKA-R4 LASER SCANNING SYSTEM
Akhka-R4DW, shown in Fig. 1, is a backpack mobile laser scanning system that can collect high precision 3D topographic data by kinematic means. The system operation is based on GNSS-IMU (Global Navigation Satellite System -Inertial Measurement Unit) positioning (Table 1). The positioning system observes GPS and GLONASS constellations for solving the instantaneous position of the platform. Precise platform movements are observed by the near-navigation grade IMU to produce attitude (Roll, Pitch, and Yaw) data for precise georeferencing of the laser scanner and image data. 3D data collection is carried out by two profiling laser scanners, each operating at different wavelengths to allow improved spectral sampling over monochromatic intensity lidar data of objects and surfaces of interest. The 3D measurement capability is achieved with two precision profiling mobile laser scanners. The primary scanner RIEGL VUX-1HA (RIEGL Laser Measurement Systems GmbH, Austria) provides for millimetre level ranging accuracy and precision. High density data is obtained using 1017kHz pulse repetition and 250 lines per second scanning rates. The scanner operates at 1550 nm wavelength giving, in addition to the 3D location of each measured point, also information of the surface reflectivity, signal amplitude and echo length. The secondary scanner, RIEGL miniVUX-1UAV operates at 905 nm. The scanning performance is moderate at 100 kHz pulse repetition rate and 100 lines per second scanning, but reasonable for backpack and other slow velocity applications. Both scanners, and the system layout, provide 360 degrees FoV (Field of View) to map the surroundings in cross-track scanning. To achieve that the secondary scanner is tilted about 30 degrees in the backpack rig to enable full FoV past the optics of the VUX-1HA (See Fig. 1A). The primary scanner is in a position for nominally vertical scanning. However, the mount on the backpack is in practice slightly tilted forwards while operating.  For image data the system is fitted with FLIR Ladybug5+ panoramic camera (FLIR, Canada), which consists of six cameras providing 5Mpix images each (Fig. 2B). Image capture is synchronized to the system to solve for camera position and orientation at the time of each exposure (cf. Fig 1B). For terrain modelling at walking speeds the image exposure was set for every two seconds (0.5Hz, an image every 2-3 meters). The image exposure is controlled by the GNSS receiver and the signal timing is recorded using the receiver's Marktime log. Image pixel size was 2 mm at 10 m distance.
Software can output each cameras' image in the sequence separately, or as a stitched panorama. The image data was not used in the terrain modelling for this paper, but remains a future effort to add RGB to models, and to study feasibility of the multi-angular imagery to produce terrain models using structure from motion (see e.g. Ullman, 1979;Johnson et al., 2014).

TERRAIN DATA COLLECTION ON AXEL HEIBERG ISLAND
The site of this study is located on Axel Heiberg Island, a desert island of 43,178 km 2 in area partially covered with glaciers belonging to Nunavut territory (79.797694N, 91.303975W). An intriguing fact is that there are petrified forests from the Eocene period on the island. The multinational field team visited numerous locations on the island during the summer season 2019 for planetary geological research and mapping. This paper focuses on the lidar system and data processing applied on data acquisition from one of those locations, specifically a thaw slump about 2 km southeast of the campsite of the expedition. Figure 2. A: The thaw slump area is seen on the opposing slope across a valley (yellow). B: The outcrop of ice exposed by the thaw slump seen by Ladybug5+ camera. The smaller slump higher on the slope is also seen as a scar to the top right.
The study site is characterized by high-centred polygons that overlie a massive ground ice deposit. A 150 m long outcrop of ice was exposed by a presently active thaw slump on a southwest sloping deposit shown in Fig. 2A, and indicated by a yellow rectangle. As described by Hibbard et al. (2020a), the polygons in the area have an average diameter of 20 m with a generally rectangular shape. The GPR readings of the soil atop the buried ground ice varies from 0.55 to 1.3 m in thickness in agreement with the estimated local active layer depth. Whitecoloured ice is visible just beneath the soil layer and is ca. 0.7 m thick overlying about 3 m of banded and fractured ice. About 1.7 m of ground ice occurs at the base of the outcrop (Fig. 2B).

GNSS-IMU trajectory processing
The trajectory required for the mapping data (lidar and images) georeferencing is typically computed in multi-pass postprocessing to provide the best possible accuracy for the solution. For that purpose a separate GNSS base station (Trimble R10 (Trimble, USA)) was erected at the field camp site (79°23'05.51 740''N, 87°51'29.18 395''W) and used to collect GPS and GLONASS satellite range observations at 5 Hz during each lidar data collection over the 11 days of the expedition. The base station was recording each day for several hours, see Table 2 for details. The position for the GNSS base station was computed using PPP (Precise Point Positioning) with precise ephemeris data.
Simultaneously on board the Akhka-R4DW, Pwrpak7 remote receiver collects the corresponding observations so that we can solve for the scanning system position as a function of time, and e.g. for errors in satellite system time and tropospheric and ionospheric disturbances and delays in the received GNSS signals during the survey. In effect we use differential GNSS to improve the positioning solution accuracy.
Fusing the IMU observations from the ISA-100C inertial unit with the GNSS data we could eventually compute precise, combined location and orientation of the mapping system. The IMU uses FOG (Fibre Optical Gyroscope) technology to observe the system rotations, angular velocities to be exact, and MEMS (Microelectromechanical System) accelerometers (along the IMU x, y and z axis) to solve for the sensor positions in high frequency supplementing the relative sparse samples of the GNSS observations. The PPP and trajectory computations were carried out using Waypoint Inertial Explorer 8.80 postprocessing software (NovAtel Inc., Canada), and tightly coupled processing. Trajectory solutions were output at 200Hz.

Point cloud processing
Processing the trajectory and raw lidar data into 3D point clouds was carried out using RIEGL proprietary RiProcess software with supplementary tools for system calibration and trajectory optimization. The data was first processed into a point cloud using initial boresight alignment angle values, and the data was filtered for long range observations and dark reflectance objects (max. range 150 m, min. reflectance -15 dB). Next in the system calibration process, the scene was automatically analysed for planar objects. Only angular parameters were solved. Translational offsets from IMU navigation centre to the scanners' origin were considered known from the design of the system with better than 1 mm accuracy. As plane search parameters for automatic tie plane search in RiProcess we used 5° minimum plane inclination, 5 cm maximum plane point deviation and required 100 points to be on each candidate plane. A 5 m search radius was used in the plane matching with 1 degree angular tolerance and 10 cm maximum plane distance. The number of observations was reduced based on the spatial distribution of the pairs and their inclination class (5° class width).
The point cloud was then transformed into a map grid projection (UTM 16) with ellipsoidal elevation (WGS84). Data was then further processed in TerraScan (TerraSolid Oy, Finland) to thin the data for elevation models.

From a point cloud to a terrain model
To turn the point cloud data into elevation model data, we needed a short pipeline. This was primarily to reduce the data to an irredundant subset of points; the original data contains in total about 6.8 billion points, 25.4GB in LAZ format. Each minute of the point cloud data from each scanner were first thinned to approximately 1 cm maximum point distance using 3D point distances for neighbourhood analysis keeping the central point within each 1 cm neighbourhood.
Next, all the thinned data were combined, but the file number ('file source ID') was preserved. For the combined data we performed a secondary thinning using a grid based method of 2 cm ground distance preserving the median elevation point within the neighbourhood. Then we filtered isolated points expecting at least 2 points within 5 cm to each other, which in general removed air points occurring due to shallow measurement angles, some long range observations, and points showing up under the ground. This provided point cloud base data for coarser resolution models. The step was performed separately for the data from each scanners.
This approach assumes that the ground is fairly bare, as was the case on the study site. In case of vegetated surfaces the processing would probably require more sophistication, and an alternative for the ground classification would be to use triangulation method, such as based on that of Axelsson (2000).

Reflectance composite
As a preliminary approach, the point cloud data was processed into reflectance composite by fusing the 3D points and lidar reflectance in to a single point cloud. The geometry of the presentation was inherited from the VUX-1HA data as is presents the most dense and accurate data. Each VUX-1HA point was associated with additional reflectance value based on the miniVUX-1UAV data by simple nearest neighbour search. However, a 5 cm distance threshold was applied to the search, and if no 905 nm data was close enough, also 1550 nm data point was discarded. Additionally, we computed NDVI (Normalized Difference Vegetation Index) for each point with reflectance from each source; R1550 and R905 (Eg. 1).
where R denotes the measured lidar reflectance. It is to be noted, that technically the 905 nm does not represent the visible red, and 1550 nm is beyond the NIR in the sense of the ordinary NDVI definition (Rouse et al., 1974), but offers a simple normalization to represent the lidar reflectance data. For colour composite the NDVI value was stored in red, R1550 in green and R905 in blue channel of the LAZ point cloud file.

PPP on the benchmark in different days
The base station observations were computed using antenna phase centre. This approach did not allow the data to be tied to existing datum and coordinate systems; that takes measurements on known points. As we did not move the base station during the campaign, for certain quality check for the PPP method we could use the solved base station coordinates for each day over the 11 day campaign to have an understanding of the absolute accuracy of the base station location.

Date 2019
Easting (  with precise ephemeris data for each of the days. The standard deviation of the daily point coordinates over the time period were 6 mm for Easting, 4 mm for Northing and 8 mm in elevation. The mean coordinates were to serve as the base station location for all the subsequent trajectory processing of the expedition data.

Trajectory computation
The base station distance to the approximate centre of the thaw slump was 1.7 km. The duration of the data collection, i.e. the Akhka-R4DW positioning data record time, was 4 hours 23 minutes, and the average speed was 3.3 km/h (0.9 m/s). During the mapping there were from 12 to 17 GPS and GLONASS satellites observable at all times. We used 15° elevation mask to filter low satellites out to improve accuracy. Fig. 3 shows the post-processing result of the trajectory along with plots of the estimated position and attitude accuracies.
The positional accuracy (see Fig. 3 for the plot) shows 5 mm 3D RMS, while the maximum 3D error was 16 mm. The RMS error for the East and North coordinates were 1.5 mm and 1.7 mm, respectively. Height RMS error was three fold, being 4.9 mm.
Estimated attitude accuracy (Fig. 3)   Combined separation is a way to analyse the quality of the positional trajectory solution. That is computed by comparing the forward and reverse GNSS-IMU solutions to reveal problems in the solution. The combined separation computed for the trajectory is plotted on top of Fig. 4. The mean separation between the solutions was 5 mm in 3D; standard deviation shows 9 mm in 3D, RMS about 11 mm. The largest deviations are in the East coordinate, which may partially be explained by the general slope orientation on the site. There is also one large peak in the separation plot, which seems to be as a result of the operator tripping over in the mud disturbing the satellite view for a brief moment of horror. The attitude separation shows some slight trend in the heading and roll angles, while the pitch angle shows steady, low deviation behaviour. The trend may be result from improper alignment in the static initialization of the IMU at the end of the data, as seen from the very last 350 s section of the data in the combined separation plot in Fig. 4. There, in contrast to the less than 2 mm separation at the beginning of the data, the plot shows about 10 mm separation for height and East directions.

Boresight calibration
In the boresight calibration we estimated the VUX-1HA and miniVUX-1UAV scanners' orientation parameters with respect to the system IMU to compute the point cloud data using ''Scan Data Adjustment' tool in the RiProcess. Robust optimization method was used, and the tie plane observations were weighted by range. We only calibrated the rotation angles as free parameters, since the translations between the IMU centre of gravity (IMU origin) and the origins of the two scanners were known from the system design at better than 1 mm accuracy.

Scanner
Roll Pitch Heading Residual 3D SD Obs.
VUX 0.07352 -0.06226 -0.17559 0.027 m 522734 miniVUX -29.82013 -0.04662 0.55427 Table 3. Boresight calibration result for the six (3+3) angular orientation parameters. Table 3 summarizes the solutions for the parameters, residual standard deviation, the number of planar matches and logarithmic residual histogram. The automated plane search, matching, and computing the calibration solution took about 2 h. The table also includes an orientation chart of the plane match pairs. The site geometry supports the calibration solution with sufficient distribution of surface orientations.
It is to be noted that in general case, the trajectory accuracy assumedly affects the boresight calibration quality. For poor positioning data, or when the scene does not contain sufficiently planar features, the solution might be weak, if not unsolvable.

Point cloud data:
The entire dataset covers an area of 1260 m by 800 m, covering also a gully to the north of the large thaw slump, Fig. 8, and a second thaw slump higher on the slope to the southeast of the one discussed here (see the trajectory in Fig. 3 for reference, and Fig. 7 for terrain detail). Based on the point cloud data, the large thaw slump was about 150 m wide, and the ice wall to the northeast, facing southwest about 10 m high. The soil overlying the buried ice outcrop ranged from 0.55 to 1.3 m thick, which is similar in range to the estimated local active layer depth (Hibbard et al., 2020a). Fig. 5 shows a top view of the thaw slump point cloud data. The slump is formed in the west-southwest inclined slope of a deposit covered with large contraction polygons. The trajectory around the slump mapped the surrounding polygon fields and the melt water channels emanating from the thaw slump excavating further the slope. The current channel is oriented towards the general slope direction with a parallel, but at the time dry, secondary flow pattern. To the northwest of the slump there was a small bond with at the time dry connection to a tertiary channel towards west. The terrain around the bond to the west shows somewhat lookalike features to so-called brainterrain type of features (see e.g. Hibbard et al., 2020b). The point density varies across the data being the densest in the polygon field around the slump itself, and close to the trajectory, about 5000 pt/m 2 . Within the slump the ground was more unevenly covered due to the sinking mud, but the bottom of the slump was also seen from the top and sides, so at sparsest, judging by the trajectory, the data shows density of about 1900 pt/m 2 . At the periphery the point density gradually drops and large data gaps become apparent due to large beam inclination and terrain undulation (see e.g. top right of Fig. 5).

Reflectance map:
The composite reflectance data in Fig. 6 show rather uniform brightness of the top soil on polygons and dried up soil atop wetter material (light brown) while mud and ground ice surfaces show significant reduction in reflectance (red tones). The fresh exposed soil on the top soil/ground ice interface show strong reflectance on 1550 nm (green). Small pooling water absorbs the laser beams completely, those are seen as holes in the data.
The top soil shows high reflectance, whereas ice, wet mud and water show low reflectance at 1550 nm wavelength. When comparing the reflectance data in the two wavelengths, the moisture has a similar dampening effect to the signals. However, for the 905 nm, the overall reflectance level is about 3 dB lower for the topsoil. On the mud plain to the top left of the Fig. 5, the 905 nm shows moderate brightening whereas 1550 nm seem to express slightly diminishing reflectance at the same locations. The banks of the flow channel do not show similar increase in the reflectance for the 905 nm in comparison to that seen on 1550 nm. The low vegetation present in the upper slump region, mainly tussocks, show low reflectance in contrast to much brighter soil. Figure 6. Reflectance composite image of the thaw slump and surrounding polygon terrain.

Elevation model:
The 3D point cloud data can be visualized using the elevation only to exhibit the topographic relations in the site. That is the approach to make analysis of the sloping, depth and width of the polygon troughs, spatial relations of features of interest, and local watersheds to name a few. With the thinning process the amount of points was reduced down to 414 million points with an approximate 2 cm point distance model, which in effect gives 2500 pt/m 2 , in practice this could be higher due to 3D shapes of the surfaces, whereas the original point cloud data was of much higher resolution. Such data gives very uniform reconstruction of the terrain features, as seen in Fig. 7. The data was also further reduced down to about 88 million points reconstructing the terrain in an irregular point distribution with approximately 5 cm point distances. Fig 8. illustrates that data as rendered into a TIN shading (Triangulated Irregular Network). The data in Fig. 9 as a traditional 2.5D DEM (Digital Elevation Model) with 10 cm grid cell (pixel) size and the elevation value expresses the minimum point elevation within each cell. Though "sparse" in comparison to the original point cloud data, the accurate and detailed sampling of the surface retains the surface characteristics suitable for many geomorphological analysis. We can, for example, detect younger, secondary polygon forming taking place on the larger polygons. The data captures also a recent mudslide in the middle of the figure. The data can also resolve the shallow ripples on the mud surfaces on the left side of the thaw slump, and retains structural fracture features on the ground ice wall itself. Figure 8. Triangulated, shaded surface model of the thaw slump area reveals micro topography of the landforms. Polygon constrained water flow patterns can also be seen in the top left of the image running towards a larger channel. Figure 9. 10 cm pixel cell minimum elevation representation of a section of the data reveals even the tiny details of the topography. An interesting feature are the small ripples across the surface indicating slow flow of material. On the polygon edges the data could reconstruct the cracking of the soil, water flow patterns, and secondary polygon forming.

DISCUSSION
The Akhka-R4DW system could be run on standard UAS LiPo batteries. This helps in having availability to resupply in any decent hobby store around the globe. 2kW gasoline aggregate suffices to recharge the batteries and tablet computer in the field. The system consumes about 116W, so it can run about 70 minutes on a single 6600 mAh battery. With low power consumption miniVUX-1UAV provides a very attractive option for lightweight, yet reasonably capable backpack scanning.
Despite the fact that in general the resolution of the ground surface data is good, the increasing beam inclination angle has its effect on data. Undulating terrain, troughs, glens, dunes, rocks and other elevation changes obstruct the view and the data suffers from elongated shadows. That is why the trajectory may seem exceedingly curvy and lengthy.
The dual-wavelength lidar data is expected to provide insights into the surface decomposition in different terrains, soil types, rocks and so on. The case study shows the both instruments provide decent and seemingly coherent reflectance data for further analysis. From the data we could make the first preliminary analysis to see what could be possible with such information in terrain modelling applications.
More detailed research is needed to characterize the scanners spectrally. It takes a method to combine the reflectance data in 3D, calibration with reflectance reference, and accounting for the range and incidence angle effects to perfect the reflectance data. But even as such the data shows great potential in characterising objects of interest.
The processing pipeline is based on tools implemented in commercial software tools. The most time consuming part of the process was to export the point cloud data from the RiProcess into LAZ-format, even using fast SSD drive, for further terrain model processing in TerraScan. There is plenty of room for improvement for RIEGL for their implementation.

CONCLUSIONS
In this paper we introduce Akhka-R4DW, the first dualwavelength kinematic backpack laser scanning system and present its application to high resolution 3D terrain modelling of permafrost landforms in the Arctic. The system is implemented on top of GNSS-IMU positioning with two profiling mobile laser scanners operating at 1550 nm and 905 nm wavelengths, supplemented with panoramic multi-camera system. The purpose of the instrument is to provide high accuracy 3D and thematic object data for environmental, industrial and planetary applications. The dual-wavelength approach is targeted for improved spectral sampling to overcome limitations of monochromatic lidar systems.
Based on the results, the proposed backpack laser scanning provides accurate and precise 3D data and performance over considerable land surface area for detailed elevation modelling and analysis of the morphology of features of interest. Using the backpack system areas covering several square kilometres could be mapped. The resulting trajectory shows millimetre level agreement; the positional accuracy shows 5 mm 3D RMS, while the maximum 3D error was 16 mm. Cumulative, trajectory related 3D RMS error estimate accounting both position and attitude uncertainties for a single lidar point at 125 m range was 13.5 mm. The residual standard deviation for the robust, observation range weighted boresight calibration of the system was 27 mm expressing the overall estimate of the point cloud geometric quality. Nonetheless, we estimate that the final elevation model data is more accurate, probably down to 8-10 mm level, as the DEM process favours short range observations due to the median elevation selection.
The accurate topography obtained allows for more elaborated and precise characterization of specific factors describing particular sites and landforms to understand and model the influence of, e.g., thermal accumulation through the general slope direction, steepness and subtle variations therein, reconstruction of surface decomposition and fracturing of the surface, troughs, faults, and channels and gullies formed by the melt water and frost, on processes in permafrost regions. The precision of the point cloud provide mm scale ground truth for interpreting of satellite radar signals from planetary surfaces.
Although not studied in this work, the obtained accuracy suggests use of the data in monitoring spatial and spectral changes over time at centimetre level in the future.
The obtained density and accuracy of the point cloud data also permits combination of the lidar reflectance data for spectral products, possibly e.g., for reconstruction of the surface composition, roughness and separation of vegetation. However, more methodological development and analysis is needed to fully explore and report the spectral behaviour on different calibration and natural surface materials and application of the data in point cloud processing, structural analysis and data classification through derivative spectral indices and ratios.