Measuring wheel ruts with close-range photogrammetry

We demonstrate the efﬁcacy of using close-range photogrammetry from a consumer grade camera as a tool in generating high-resolution, three-dimensional coloured point clouds for detailed analysis or monitoring of wheel ruts. Ground-based timber harvesting results in vehicle trafﬁc on 12–70 per cent of the site, depending on the system used, with a variable probability of causing detrimental soil disturbance depending on climatic, hydrological and soil conditions at the time of harvest. Applying the technique described in this article can reduce the workload associated with the conventional manual measurement of wheel ruts, while providing a greatly enhanced source of information that can be used in analysing both physical and biological impact, or stored in a repository for later operation management or monitoring. Approaches for deriving and quantifying properties such as rut depths and soil displacement volumes are also presented. In evaluating the potential for widespread adoption of the method among forest or environmental managers, the study also presents the work-ﬂowandprovidesacomparisonoftheeaseofuseandqualityoftheresultsobtainedfromonecommercialandtwoopensourceimageprocessingsoftwarepackages.Resultsfromacasestudyshowednosigniﬁcantdifferencebetweenpackagesonpointcloudqualityintermsofmodeldistortion.Comparisonofphotogrammetricproﬁlesagainstproﬁlesmeasuredmanuallyresultedinrootmeansquareerrorsofbetween2.07and3.84cmforﬁveselectedroadproﬁles.Maximalwheelrutdepthforthreedifferentmodelswere1.15,0.99and1.01m,andesti-matedrutvolumeswere9.84,9.10and9.09m 3 , respectively, for 22.5 m long sections.


Introduction
Timber harvesting is an expensive part of forest management and has the largest potential for generating untoward environmental, ecological and aesthetic impacts (Duncker et al., 2012). Mechanized ground-based harvesting systems account for the vast majority of timber harvested commercially across the world, with the risk of soil substrate damage in the form of compaction, shearing or rutting being greatest during the extraction phase. Extraction involves the use of wheeled or tracked machinery ranging from agricultural tractors through bulldozers, skidders and forwarders, which either drag or carry timber from within the stand to the forest road network. Forest certification standards and the implementation of sustainable forest management practices advocate the avoidance or amelioration of wheel rutting from forest harvesting operations.
Wheel rutting occurs when the ground pressure from the tire or track exceeds the bearing capacity of the soil, causing a compaction rut which, under increasing tractive demand, soon develops into a shearing rut due to wheel slippage (Eliasson and Wä sterlund, 2007). Many technical developments and management considerations have been given to mitigating the impacts of forest operations on soils (Talbot and Astrup, 2014). However, the future effects of climate change, including longer periods of saturated, unfrozen soils in the boreal zone (Soja et al., 2007), imply that wheel rutting, which is already a widespread problem (Cambi et al., 2015), is likely to become a more critical issue in the future. Even in the best managed operations with sophisticated harvester/forwarder cut-to-length systems, vehicle traffic occurs on at least 12 per cent of the site (Eliasson, 2005).
Wheel rutting has several documented negative environmental effects. It has been shown to be detrimental to tree growth due to both physical root damage and reduced soil porosity (Wä sterlund, 1994;Quesnel and Curran, 2000;Sirén et al., 2013). A continuous linear wheel track leads to waterlogging across the slope or channelling if perpendicular to the slope (Startsev and McNabb, 2009). Rutting and other soil disturbance from machine traffic can also promote transmission of pathogenic fungi through root networks (Thor and Stenlid, 2005) and mobilize heavy metals, although the direct cause and effect has not been fully quantified (de Wit et al., 2014;Eklö f et al., 2014). Rutting can also exacerbate oxidation and the emissions of greenhouse gases from otherwise stable forest soils (Ampoorter et al., 2010) and ruts can impede passage to recreational users of the forest and are widely considered to have a strongly negative aesthetic impact (Gundersen and Frivold, 2008).
Wheel ruts have traditionally been manually quantified for both research and management purposes, where the horizontal and vertical displacement is determined using a levelled hurdle and ruler, measuring tape or infrared range finder (Nugent et al., 2003). The labour demanding nature of such measurements implies that only a relatively small number of transects can be sampled, making it a task well suited for the application of modern geospatial technologies (Koren et al., 2014).
Work has been done on ruts and deflection measurements with vehicle mounted laser rut measurement systems on paved public roads (Laurent et al., 2012). In such applications, the high cost of these systems is offset by the speed of assessment and socioeconomic importance of their application to public infrastructure, neither of these criteria being valid for dispersed forest harvesting sites. Similar research has also investigated using airborne laser scanning (ALS) in order to address consequences of land use on soils, specifically the mapping of gullies and headwater streams (James et al., 2007;Perroy et al., 2010). ALS data from regional, multipurpose scans contain too few points (2 -5 pts m 22 ), and the data are captured too infrequently at stand level (typically only once), for the modelling of sporadically occurring structures such as wheel ruts, while higher resolution drone-based laser scanning is only becoming available as a result of continuous reductions in the size, mass and data storage demands of laser scanners. Koren et al. (2014) used terrestrial laser scanning from a static set-up in accurately assessing rutting depths and surface disturbance after harvesting with a high-density point cloud ( 50 points m 22 ).
Other methods under development for assessing wheel rutting during operation include the monitoring of the machine power requirement from the CAN-bus data as a proxy in estimating rolling resistance on the machine, corroborated by simultaneous ultrasonic distance ranging (Ala-Ilomä ki et al., 2012). Another method estimates wheel slip, which can also indicate rut development by comparing CAN-bus information on wheel movement with real-time kinetic (RTK) positioning (Ringdahl et al., 2012). Unfortunately RTK GPS is expensive to install on forest machines and there are still positional accuracy issues to resolve.
Photogrammetry is a method of remote sensing that offers a solution to many of the above issues and has been applied to soil surface modelling and pavement distortion quantification in various forms for some time (Warner, 1995;Obiadat et al., 1997). The use of close-range photogrammetry in mapping soil surface structure was demonstrated more than 20 years ago (Warner, 1995), but developments in consumer grade cameras and improved digital workflows have simplified the derivation of high-resolution surface models (Turner et al., 2012) and have been used in applied contexts in forestry (e.g. Pierzchała et al., 2014a, b) and agriculture (Nouwakpo and Huang, 2012). An advantage of photogrammetry-based methods using modern compact cameras is that they can be used handheld or mounted on forest machines or low-cost drones, while laser scanners remain comparatively heavy and costly. Close-range photogrammetry provides very high-resolution point clouds. Efficient threedimensional (3D) reconstruction methods from images have been made available in numerous web-based processing services, which provide good quality textured 3D models that can be exported. However, these web-based processing services do not allow for additional user input such as intrinsic and extrinsic camera parameters, and therefore often produce somewhat distorted scale models in arbitrary coordinate systems. These models are primarily designed to create sufficiently accurate 3D representation of the objects and are typically not intended for retrieving accurate spatial information.
Given the importance of reducing wheel rutting for sustainable forest management, the development of efficient methods for measuring, describing and disseminating information on their current prevalence and severity is an important step in mitigating anticipated rutting damage in the future.
The aim of this article is to present a method for remotely retrieving spatial information on wheel rutting that provides a low cost, precise and accurate estimation of the rut by means of close-range photogrammetry, 3D modelling and the analysis of meshed surfaces. A secondary aim is to compare the quality of results from two open source web-based processing services with commercial modelling software in order to evaluate the potential for a more widespread adoption of the method in forest management.

Field trial
The trial was carried out on a site on the eastern shore of Oslo fjord, Norway (Lat. 59.611558; Long. 10.697114) in the late summer of 2014. A clearcut of 1.5 ha had been carried out and the timber extracted in the spring of 2014. The stand was dominated by Norway spruce (Picea abies) and growing on a flat to slightly undulating, partially waterlogged moraine with intermittently exposed bedrock. The machine used for timber extraction was a Komatsu 860.4 forwarder fitted with steel band tracks on the bogie wheel sets. While wheel rutting was apparent on a large extent of the site, only a portion of the exit/entrance to the site of roughly 22.5 m in length was selected for the analysis on the basis of this area showing severe rutting, its easy accessibility and sufficient extent in testing the method ( Figure 1). Due to the length of Forestry time elapsed since the harvest (one summer), shrubs, grass and harvesting residues needed to be cut and removed in a corridor of 5 m in width on either side of the wheel ruts. This was necessary in providing a clearer view of the original terrain surface.
Close-range photogrammetry data were collected using a consumer grade compact camera (Panasonic DMC-SZ8) suspended from a manually carried 4 m long pole (Table 1). Image acquisition was done using a Wi-Fi camera control App installed on an iPad, which allowed the operator to view and control for sufficient image overlap and quality. The area covered in the trial was 210 m 2 . Five checkerboard ground control points (GCPs) were positioned on each side of the ruts, and at a lateral distance far enough from the ruts to ensure that the terrain surface had not been deformed in any way by the forwarder passes. Edge and diagonal distances between subsequent GCPs were measured with a steel tape and used in scaling the reconstructed surfaces ( Figure 2). Once the data were processed into the final textured mesh, the GCPs were visible in the textured mesh as well as at singular points in the point cloud.

Surface reconstruction
A total number of 56 images were processed using three different software packages: Agisoft PhotoScan w , and the online web services cmpmvs w and Autodesk 123D Catch w . These software packages apply a Structure-From-Motion (SFM) process that relies on feature-based image matching in order to estimate camera poses and thus retrieve full 3D scene reconstruction (Wu, 2011).
Agisoft PhotoScan w is commercial software that allows for detailed model parameterization (Agisoft, 2012). It takes a series of images as an input and estimates camera intrinsic parameters. The software facilitates direct or indirect geo-referencing either by introducing the positions of cameras and rotation data or by GCPs. Furthermore, it allows for multiple exports of RGB point clouds with normals, textured meshes and digital terrain models in raster format. The camera was pre-calibrated using the Matlab w camera calibration toolbox. cmpmvs w is a software developed at the Centre for Machine Perception, FEE, CTU, in Prague. It has an image processing SFM Web Service and MVS -Multi-View Reconstruction Software (Jancosek and Pajdla, 2011). Image matching is supported by a Bundler tool (Snavely et al., 2006) or VisualSFM (Changchang, 2013).
Autodesk 123D Catch is one of various web-based services for image processing. As an input, it requires a series of images and returns a textured mesh. Autodesk 123D Catch is used without prior calibration and relies solely on self-calibration procedures. Autodesk does not provide the possibility of introducing a camera calibration matrix prior to the processing. All the created point clouds and meshes were further processed with the aid of Cloud Compare w , Meshlab w and Matlab w .

Scaling and quality measurement
After processing the point clouds with the three software packages (Photo-Scan, cmpmvs and 123D Catch), the three resulting surface models had to be scaled to convert them from arbitrary dimensions and allow for measurement of true distances. The surface models were adjusted using uniform scaling based on the distance measurements between the GCPs (10 nodes in Figure 2) and linear distances (13 measurements shown in Figure 3). The corresponding points in the point cloud were selected by identifying the position of the GCPs on the textured meshes ( Figure 3). The mean scale factor was calculated as: where S is a mean scaling factor, d i is a real (measured) distance and d ′ i is the point cloud distance. Mean scale factor is then used to standardize individual scale factors for evaluating the deviation of each observation. This is done to verify whether the output surface models preserve distances and do not exhibit significant model distortion. It also shows that manual point selection is sufficient in scaling the model when working in an arbitrary coordinate frame.
Following this, a similarity transformation was carried out with 7 parameters using the 10 control points that recorded DGPS positions, transforming all the surface models to the same coordinate frame. The targeted coordinate system that was used is altered ETRS89/UTM zone 32N (EPSG: 25832) where the first four digits on Northing and three digits on Easting were truncated. This was done because mesh processing software performs better with shorter digits, and the practice is recommended in avoiding redundancies. Truncating also significantly reduces the file size and thus speeds up the processing time.
Model quality was assessed separately for the terrain models generated by each of three software solutions. Quality was defined according to the surface model distortion caused by insufficient self-calibration as it is known that methods using manual camera calibration outperform selfcalibration techniques (Mendonça and Cipolla, 1999). This distortion is mainly radial and results in model deformation known as the 'bowl effect'. By assuming marginal geometrical distortions, an Iterative Closest Point (ICP) algorithm was run to finely register the models. ICP aims at minimizing the alignment error (E) with respect to the rotation (R) and translation t: where (p i ,q i ) are collections of points with normals n i . The reference model used here was the surface model from precalibrated cameras using PhotoScan. After this, the signed distance from each point to a reference mesh was computed, resulting in a point cloud with assigned new scalar values representing the signed distance. Signed distance is essentially a distance between a point and a plane. Given a plane: and a point x 0 ¼ (x 0 ,y 0 ,z 0 ), the normal vector to the plane is given by:  Figure 2 In-field measurement of distances between GCPs (filled circles) used in scaling the images. The numbers indicate the sequence in which the measurements were done.
Measuring wheel ruts with close-range photogrammetry And a vector from the plane to the point is given by: Projecting w onto v gives the distance D from the point to the plane as: Dropping the absolute value gives the signed distance: which is positive if x 0 is on the same side of the plane as the normal vector v and negative if it is on the opposite side (Weisstein, 2015). All three models were clipped to one common bounding box, and subsampled using a random subsampling method with a fixed target with 10 000 points. The distribution of signed distance was then studied.

Comparison with conventional registration
Wheel rut profiles have traditionally been measured in transects using a horizontal hurdle and measuring rod. The same method was used as a control in this study, where five transects of between 5 and 6 m were measured manually at centimetre vertical accuracy and at 25-cm intervals across the transect, and these points were assumed to represent the true profile ( Figure 5). The photogrammetry-derived surface models were then compared against these manually measured profiles by registering corresponding point sets from both origins. Initially, 30 cm wide strips were clipped from an aerial view of the generated terrain model for every transect then rotated into cross-sectional profile position using four points (two GCP and two from vertical manual measurements). After transformation, the Forestry Z values (depth values) were removed in order to reduce the problem to two dimensions. The resulting point sets were then fitted against the point set from the manual measurements. The ICP algorithm was run to register two point sets by minimizing the root mean square (RMS), with 15 iterations. The ICP returned the transformation matrix that was then applied to the manual point set.

Depth and volume calculations
In order to calculate the volume of the displaced soil, the original terrain surface needed to be reconstructed. The surface that had been altered by the wheel rutting was visually identified and manually clipped from the surface model. The remaining part was used for surface reconstruction. For our assessment, we chose Poisson Surface reconstruction (Kazhdan et al., 2006), implementing adapted octree (Meagher, 1980). In this case, octree level 28 was used as it produces a fine resolution result. Only the PhotoScan-derived model was used for computing the reconstructed (prerutting) surface. The extraction trail model clipped from the larger surface model was compared with the reconstructed terrain surface (pre-rutting) and the signed distance was computed and assigned as a scalar field to the corresponding vertices of the modelled surface. Positive values were assigned to all vertices that were located above the reconstructed surface and negative values are assigned to all values under the reconstructed surface.
The next step was to filter out the negative values by thresholding the signed distance values. The filtered mesh could then be fused together with the reconstructed mesh, forming a watertight 3D surface. After fusion, we checked whether the resulting mesh was manifold and then the volume calculation was performed. The manifold check is easily carried out in Meshlab w . The volume was calculated using a polygon mesh voxelization that uses the ray intersection method of voxelization (Patil and Ravi, 2005). This outputs a 3D logical matrix with n-elements defined by the user, which determines the level of accuracy required. As input variables, we provide the number of elements for every dimension in the matrix; however, voxels will not be cubic but cuboid. In order to calculate the total soil volume, we apply the following equation: where v is a voxel, b is a three-element vector representing the dimensions of the bounding box for the input mesh, g is a three-element vector representing the number of elements in logical matrix and n 1 is count of true values in the logical matrix. Volume results are reported in m 3 . Voxelization was done for three different meshes with a reference-reconstructed mesh from the PhotoScan model.

Results
The image processing resulted in three 3D models of surface that contain 3D point locations with normal, meshed surfaces and texture (Figure 3). The output models differ slightly in terms of their relative size. Surface texture is well modelled, and the continuous caterpillar track marks are easily distinguishable.

Scaling
After scaling the point clouds, we obtained different scaling factors as an average of derived scaling factors for all 13 measurements in each point cloud. After applying an averaged scaling factor for all measurements, we retrieved the real dimension by dividing by an original distance (Figure 4).
The original average scale factors for the three models were as follows: PhotoScan ¼ 48.22;123D Catch ¼ 50.99;cmpmvs ¼ 83.22. Figure 4 shows similar trends in terms of the sign of the standardized scaling factor both for uniform scalar multiplication and similarity transformation. This could be explained by model distortion that reaches maximal deviation at the end of the model. Big positive values correspond with measurement numbers 1, 8 and 9 that are located at the beginning of the section and big negative values 3, 4, 5 and 13 are located at the other end (compare with Figure 2). An alignment comparison of two meshes against the reference surface shows only small standard deviations (Table 2).
A comparison of the horizontal hurdle and measuring rod manually measured profile with the photogrammetric records showed good correspondence. Comparison of the manual measurements and the PhotoScan model resulted in root mean square error (RMSE) values for each of the five transects were 1 -2.73, 2 -2.71, 3 -2.07, 4-3.84 and 5-3.63 cm, respectively ( Figure 5). This value is satisfactory taking into account the manual measurement method precision.
The reconstructed model clearly showed the upwelling and rutting which resulted in overall soil displacement ranging from positive values of 60 cm (above original surface) to roughly 1 m in depth ( Figure 6).
The three software packages gave similar though not identical density distributions in representing rut depth, where the modus Figure 4 Standardized scaling factors for each of the 13 field-based distance measurements between the GCPs. Values for PhotoScan, cloud123D catch and cloud cmpmvs stem from uniform scaling based on linear distances between GCPs. The measurement numbers (x-axis) correspond to the distance measurement numbers illustrated in Figure 2. Measuring wheel ruts with close-range photogrammetry for the PhotoScan model lay slightly higher than for the two others (Figure 7), and are summarized in Table 3. They present a maximum rut depth of 1 m (min.) and soil accumulation of roughly 50 cm (max.). Volume calculation returns a 3D matrix with 50 × 50× 50 dimensions and total voxel count of 125 000 (Table 4). The bounding box is automatically derived from input STereoLithography mesh. Voxel dimensions are not cubic and their size explicitly relates to bounding box size. Voxel count describes how many voxels fell into the model. Full voxel fraction shows the share of full voxels (representing the volume of ruts) to a total number of voxels felling into the bounding box.

Discussion
Close-range photogrammetry provided a method for quantifying the dimensions of wheel ruts after timber harvesting that is at least as accurate as the current manual method of measurement and considerably less labour intensive. The comparison of the horizontal hurdle and measuring rod manually measured profile with the photogrammetric records showed good correspondence with low associated RMSEs. Part of these errors must be expected to arise from manual measurement inaccuracies from, for example, placing the measuring rod on a stone or in a small hollow. Nouwakpo and Huang (2012) compared photogrammetry results of rut profiles with high-density terrestrial laser scanning and found RMS errors related to height (Z) of only 6 -8 mm, which is within the same range as the results presented here. The rut is described as a continuous hull, not only in linear transects, and therefore enables both length and volume to be calculated accurately. Establishing accurate and specific data on volumes can be useful if they are to be monitored for change (Oleire-Oltmanns et al., 2012).
The presence of vegetation or harvesting slash can cause problems in estimating true rut depth as photogrammetry provides only a surface and not a terrain model. The long duration (one summer) between harvesting and measurement in our case meant that grass and shrubs had established themselves adjacent to the tracks and had to be mowed in order to get a good representation of the terrain surface. However, if this method was to be applied immediately after the forest operations, issues related to vegetation should be minimal and no vegetation control would be required. The colourized point cloud record of the wheel rut that photogrammetry provides is intuitively understood and can be geo-referenced and used also for other purposes, such as vegetation monitoring (Saarinen et al., 2013).

Scaling
Scaling can be the most work demanding requirement of using this approach in the field. Through an evaluation of scaling efficacy, we could assess the limits of minimum scaling information necessary in obtaining real distances. Differences in the error distribution were caused by several factors. Inaccurate measurements of distances in the field could have been caused by a measuring tape that did not sustain linearity when suspended over longer distances. Secondly, inaccurate selection of pixels corresponding to GCPs in the images can be a source of inaccuracy. The considerable disparity observed for corresponding observations might be explained by model distortion (Figure 4). The average was low in most cases except for the cmpmvs point cloud, which reached higher levels of error due to more significant distortion (Figure 4). However, we Forestry consider the observed error levels to be acceptable for our application, which is traditionally characterized by a large measurement error and mostly measured for management purposes where millimetre accuracy would not be required.
For DEM extraction, a mandatory procedure is the acquisition of control points, which involves the use of survey equipment such as total station or RTK GNSS. However, the literature shows alternate solutions for reducing the workload. One method is to acquire X,Y,Z positions of control points using known information about the object space and the calibration frame to compute the relative position of the control points (Nouwakpo and Huang, 2012). This method uses the algorithm of Van Den Heuvel (1997) that uses parallel lines in the scene to retrieve the exterior orientation of the camera.
Another solution that can be found to minimize the workload is to obtain model dimensions using calibrated extrinsic camera parameters using wide baseline stereo vision (Antonello et al., 2013). In that method, it is shown that wide baseline stereo can significantly reduce the depth estimation errors and can result in good quality 3D dense reconstruction using, for example, Semi Global Block Matching algorithm (Hirschmü ller, 2008).
Since our method predicts wheel rut depths in arbitrary coordinate system, it is sufficient to utilize extrinsic parameters obtained from stereo rig baseline distance to obtain real distances.
The proposed workflow can be tested on a greater scale, but precision issues arise. Environmental evaluation of a large harvesting site would require the use of an UAV platform with a low elevation flight and high image resolution, which in turn produces very big datasets that are difficult to handle.

Mesh registration and signed distance distribution
The mesh registration turned out to be very accurate, and improved the initial alignment significantly. Distribution of signed distance is normal and the values of mean and median are within 1 cm (Table 2). However, some points exhibit high values of up to 15 cm (Figure 7). These are likely points with incorrect estimates of 3D position and can be treated as outliers. Some points are simply wrongly estimated due to poor textures, in this case due to surface water. Water would be commonly encountered in studies of wheel ruts in forestry, but only becomes an issue if measures on the absolute depth and not just the existence of the rut beyond a certain threshold are required. PhotoScan generated points from surface water images, but this resulted in an overestimation of their position that led to incorrect results. The problem of water in the rut could be partially mitigated by considering seasonal fluctuations in precipitation when planning to image the sites.
Signed distance is a measure that gives a direct representation of wheel rut depths. It is easy to generate when having a reference   Measuring wheel ruts with close-range photogrammetry mesh and produces easy to interpret results. Result of signed distance analysis towards the reconstructed surface (Table 3) returns similar values for three meshes except a minimum value that is 15 cm greater in the case of PhotoScan, here due to the effect of incorrect estimates of water surfaces as discussed before. Finally, we evaluate the mesh voxelization method and its accuracy by computing soil volume for three different meshes. All three meshes have similar bounding box size and are equally clipped. It turned out that the volume estimate was very similar for the three models, and is also precise, taking into account the coarse size of the voxels. The resolution of the voxel grid could have been increased, although the default setting was considered satisfactory. Differences can be explained once again by incorrect point estimates resulting in 'deeper' ruts for the PhotoScan model.

Conclusion
The outlined approach can precisely quantify the geometry of soil disturbance, but free water in the rut is an important factor that may influence the accurate implementation of the method. It is important that poorly textured surfaces are excluded from the analysis so that no incorrect estimates are included, implying that it would not be possible to get depth information under water.
Surface reconstruction methods are very accurate, but more attention must be put into reconstruction processes and validation of the results. Ground vegetation is a limitation in the reconstruction process due to occlusion of the real terrain surface, but this should be a negligible issue if the approach is applied within a short time from the forest operations.
For implementation of the method, the workload must be reduced by getting the camera translation for use in scaling. This could be done using a calibrated stereo rig, where camera translation from each stereo pair is known. By using the SFM workflow on images where image pairs with a known baseline can be selected, we can apply scaling from estimated camera poses. However, commercial SFM web services do not provide camera extrinsic parameters. A stereo rig could be used for a different type of process, namely creating depth maps from calibrated stereo views and evaluating profiles to get the information about disturbance.
The approach outlined in this study is easy to perform with the use of open source resources and consumer grade cameras and can, therefore, be applied for small-scale evaluations or as parts of sampling schemes for evaluating the extend and severity of wheel rutting.