ABSOLUTE COREGISTRATION AND DOMING CORRECTION IN ADVERSE CONDITIONS, OR HOW TO RETRIEVE SNOW DEPTH FROM DRONE FLIGHTS

Abstract. Doming, a superimposition of a scene-wide dome-shaped bias to the elevation data, is a classical problem in photogrammetry, especially when using the auto-calibration and bundle adjustment typical of the Structure-from-Motion approach. A number of solutions have been proposed, but most rely on a certain degree of freedom in the data acquisition strategy, such as the ability to acquire many ground control points or to acquire additional imagery over the area of interest with different flight parameters.In this paper, we present a method to solve for the dome when working on a snowfield in the Norwegian mountains. That specific situation leads to difficulties to fly the drone optimally and acquire ground control points due to short weather windows, high winds, cold temperatures, low sun angles and short day duration (i.e. high latitudes). Here, we use natural landmarks and pre-existing knowledge about the bare ground topography for removing the inevitable dome.



INTRODUCTION
The use of photogrammetry in geosciences has been booming in the last decade through the use of structure-from-motion (SfM), especially since the advent of relatively cheap, high quality, reliable Unmanned Aerial Vehicles (UAVs) (Koenderink, Van Doorn, 1991, Westoby et al., 2012, Boesch et al., 2016. The presence of systematic broad-scale error like doming is however still posing issues to many users. The term doming refers to the presence of a bias in the elevation superimposing a dome to the elevation data after georeferencing. It is caused by a wrong estimation of the camera calibration (James, Robson, 2014) leading to erroneous relative external orientations. In particular, the auto-calibration and bundle adjustment processes can wrongly estimate the focal length, and compensate with the radial components of the distortion model. The tie points reprojection errors can therefor be very small while the relative orientations are affected by doming. A shorter (respectively longer) calibrated focal than reality creates a dome with the center of the scene higher (respectively * Corresponding author lower) in elevation than reality, and the edges lower (respectively higher). Fig. 1 shows an elevation profile representation of this effect.
Multiple approaches have been proposed to solve the issue, for instance : • using more geolocation data to re-estimate the calibration with either more Ground Control Points (GCPs) than the strict minimum, or the inclusion of high-quality in-flight geolocation systems; • modifying traditional flight planning to counteract the issue (flying with the camera off-nadir, combining flights at different heights, adding non-nadir concentric groups of images to the bundle, etc.) (James, Robson, 2014); • increasing the complexity of the distortion models (Tournadre et al., 2015).
These options might, however, be impractical or not perform well enough for the user's needs, as shown in our case study where UAV flights are used for mapping snow depth.

MOTIVATIONS AND CHALLENGES
The study of snow cover and snow distribution is crucial both to the understanding of water resources and to the safety of people and infrastructures in mountainous areas. Frequent high resolution snow maps are important observation to validate or assimilate to snow modeling studies but keeping data acquisition reasonably cheap offers unusual challenges. The recent development of relatively low cost UAVs sprung interest by the snow scientific community to use them for photogrammetric surveys of snowfields (Koenderink, Van Doorn, 1991, Westoby et al., 2012, Boesch et al., 2016. Performing a state of the art photogrammetric survey has a number of requirements such as weather conditions (light condition, wind, temperature) and reliable and ground control that are hardly compatible with high frequency repeat survey of snow fields in mountainous natural reserves. Therefore a slew of challenges are encountered when acquiring data for snow mapping at high latitudes.
The first is that snowfields are often close to be homogeneously white, lacking the necessary texture used by the photogrammetric process. The success of snow depth studies using digital photogrammetry techniques has shown that, given the right illumination conditions, it is possible to get around this issue (Nolan et al., 2015). This, however, still requires to fly with clear skies so that the snow micro topography can be cast into shadows, providing texture to the images. Of course, this also requires overlapping images to be acquired shortly after each other to avoid shadows (or the snow itself) moving significantly between overlapping lines of flight, hence limiting complicated flight paths or very long lines of flight. Fig.2 shows what snow structures look like in sunlit and shaded areas. A is an area that is shaded by the relief to the top right, and therefore has absolutely no texture, while B is well lit from an angle and the micro-topography is highlighted. The C areas are example of stable terrain that protrude over the snow, and can be used for the corrections presented hereafter.
The second is that winter days are short at high latitudes, with less than 6h of daylight at 60 • N, and very low sun angles, leading to little light available and reduced time window to fly a drone.
The third is the rapid weather variations, an additional constrain on the already short time window available. Light drones cannot fly in high wind, and have shortened battery lifetime in cold temperature, requiring multiple flights for even relatively small areas (e.g. < 1 km 2 ). The installation of GCPs in rapidly changing weather becomes difficult, with practically no chance to prepare the site in advance. Markers are indeed at risk of being blown away by wind or being covered by snow with winds as faint as 5 m.s −1 . For this second reason, proper GCPs are rarely used in these surveys.
The fourth challenge, specific to the location, is the mobility around the site of interest within safe conditions. Based on the mean of transportation, the snowpack conditions (e.g. ice glaze, avalanche risk), and the preparedness of the operator, access to the whole area might be limited in space and time to deploy and retrieve GCPs spread homogeneously over the entire survey area.
As an additional challenge, but also because of the opportunities they offer, we only consider light quadcopter UAVs as acquisition platforms. They are easy to operate, easy to transport through rugged terrain, and have a lighter legal framework to operate (depending on country's regulations). Finding workarounds when data acquisition can't be done in a state-of-the-art manner is a common activity in low-cost surveying, as demonstrated by (Girod et al., 2017) or (Whitehead et al., 2013). Here, we present an iterative approach that attempts to use the maximum of information from all available data. In our case, the data available need to be:

OUR METHOD
• Images from the drone; • Rough georeferencing of the drone imagery from the onboard code based GNSS receiver; • A summer, snow free digital elevation model (DEM) and the associated orthoimage, both at the same ground sampling distance (GSD) as the expected product from the images.
The process, shown in the diagram in Fig.3 starts with a standard, automatic structure-from-motion photogrammetry (SfM) processing approach, finding tie points, computing a relative orientation and auto-calibration, using the GNSS data to georeference the model and perform a last bundle adjustment and a re-estimate of the calibration before the dense correlation, orthorectification and mosaicking.
Then, the external DEM is used to pick GCPs. In this first iteration, the points picked are easily recognizable, like large erratic boulders, or stable man made objects if any are visible in the imagery. The position of each point is recorded from both the external data (XY Z external−iter1 ) and from the initial DEM and orthoimage (XY Zinit). Then, using the orientation data, the XY Zinit points are back-projected in the images (IJiter1). The XY Z external−iter1 -IJiter1 pairs are then used as GCPs to refine the absolute orientation. A second set of DEM and orthoimage is then produced.
Because the GCPs are selected on protruding objects during the first iteration, the Z coordinates can be affected by the DEM regularization, or that a small error in the XY coordinate results in a relatively large error in the Z coordinate extracted from the DEMs. Therefore, the exactitude of the re-projected image coordinates is potentially affected. For this reason, a second, very similar iteration, is then performed, selecting random points (XY Z external−iter2 ) on stable ground where the terrain is softer (low absolute curvature), to reduce the error in the Z estimate. These points do not need to be particularly identifiable as the XY registration resulted from the previous iteration (iteration 1). Therefore, selecting snow-free (stable) ground is sufficient. The XY Z external−iter2 -IJiter2 pairs are used like in the first iteration, and an optimised DEM and orthoimage are then produced.
The output of the second iteration can be considered as the best product that can be made by forcing the orientation to an external product, but can often still be imperfect, especially if the external data is itself no free of doming. Because of the geometry, if the external data is still domed, the new data can't be fitted perfectly with a forcing on the orientation, and a more rudimentary approach is necessary: computing a polynomial solution to the remaining combined dome. Typically, a second order polynomial is chosen, see equation (1). The GCPs used in the second iteration (XY Z external−iter2 ) can be used a second time to compute the correction surface (see Fig.4 for the example in our use case), combined with their paired values in the DEM from the second iteration (XY Ziter3). The resulting DEM is mostly free of coregistration error with the external DEM and can therefore be used to compute elevation changes, which in our case is equivalent to snow depth.
where a − f = polynomial factors X, Y, Z = point coordinates In short, the method combines an automatic SfM routine (see INITIALISATION in Fig.3), an attempt at solving the doming analytically (see ITERATION 1 and ITERATION 2 in Fig.3), and adds a final forcing to the elevation model, removing the remaining dome (see ITERATION 3 in Fig.3).

USE CASE
The method was developed to improve data quality of snow mapping surveys conducted over a square kilometer of mountainous terrain close to the Finse Research Station in central Norway (see Fig.5) sitting at about 1200 m above sea level. As snow covers the ground for up to 7 months a year, the area of Finse offers an ideal site to study mountain hydrology, and in particular snow redistribution processes by wind. Going from an even snowpack to an uneven snowpack (shallow snow adjacent to deep snow) has profound impacts on hydrology, ecology and climate, driving the scientific need to better understand the spatial distribution of snow.
For our use case, we acquired 300 pictures with a DJI Mavic Pro quadcopter UAV flying at approximately 120m over the ground on March 3rd 2019. From the images we generated an initial DEM with a 20cm GSD. The external DEM and orthoimage were acquired in a near snow-free period with a Nikon 1 on-board an octocopter on September 21st 2015 (also at 20cm GSD). Additionally, we dragged on the snow surface by skis a sled equipped of a differential GPS receiver on the same day as the image acquisition. This dGPS profile is an independent ground truth of snow depth for evaluating snow maps produce by the method presented herein. The 2015 and 2019 orthoimages as well as the ski track are shown in Fig.6. The area exhibits a series of ridges and troughs where snow is either blown away or accumulated depending on local wind speed (Fig. 2). We therefore used the ridges to control and estimate the magnitude of the dome in between the iteration steps.
The initial processing resulted in a DEM that is affected by a strong dome (see Fig.7A), as well as a ca. 65m vertical shift compared to the reference DEM because of a mismatch in the vertical reference (Ellipsoid vs Geoid). The first and second iterations of GCP picking on the reference DEM (see the GCPs in Fig.6) led to a better fit in all dimensions, but failed at correcting the dome to an acceptable level (see Fig.7B).
Finally, fitting a second degree polynomial in XY (equation (1)) corrected the remaining dome. Fig.4 shows the fitted surface and the points used to estimate the best fit via a simple linear regression. Over the area of interest, the fitted polynomial correction surface reaches extremes as low as −1.7m and as high as +2.5m. Given that the signal sought is of a similar order of magnitude, and the local noise is in the order of 0.1m, correcting for the dome has a significant impact.
From our final product, we found snow depths ranging from 0 to 6 m over the entire area (see Fig.7C), to the exception where snow was present at the acquisition of the near-snow-free DEM. In those areas, generation of a DEM failed in the near-snow-free acquisition since the images' exposures were saturated on snow surfaces.
We compared our values to an in-situ dGPS track obtained by dragging a receiver over the snow surface behind a snow-scooter. The final overall linear regression between the snow depths from the corrected DEM (iteration 3) and the dGPS track points yield a r 2 of 0.83 (Fig.7D). With the snow depth derived from DEM produced at iteration 1, this same r 2 drops to 0.47. So iteration 2 and 3 improved significantly the quality of snow depth measurements from this simplistic winter drone survey. The errors are not all coming from the UAV data, but also from the dGPS, as the terrain is very rugged. A small error in horizontal positioning translates to a bias in the estimated snow depth.

IMPLEMENTATION
The method was tested and developed using scripted calls to the free open source MicMac photogrammetric library (hosted at the French National Institute of Geographic and Forest Information -IGN - (Rupnik et al., 2017, MicMac Development Team, 2020). The iterative process described in Fig.3 is available at Figure 6. Map of the snow field used for this study with the picked GCPs and the validation dGNSS track. (Girod, Filhol, 2020). The final step of the process (polynomial fitting of the remaining dome) was implemented both Python and in MicMac (as mm3d PostProc Banana), named after one of the colloquialisms of the doming effect.

CONCLUSION
This study shows that given imagery with sufficient texture, an external snow-free DEM of reference, and stable landmarks, it is possible to obtain a high resolution snow-map product at a high GSD, even with a low-cost UAV and no time-concurrent ground control. The reference DEM doesn't even need to be free of internal distortion either, as the computed dome correction fits the combine doming, which is sufficient when the desired final product is the elevation difference between the two models. While the final step of the process is admittedly a rather crude polynomial fit, it does improve significantly the quality of the final product.

ACKNOWLEDGMENTS
This study was funded by the University of Oslo eInfrastructure Competence Hub Geohive as well as the department of Geosciences of the University of Oslo.