Spatial Co-Registration of Ultra-High Resolution Visible, Multispectral and Thermal Images Acquired with a Micro-UAV over Antarctic Moss Beds

2.1. Test Sites

In Turner et al. [26] we introduced a technique to georeference and mosaic multiple visible images collected by a MUAV. In the present study we are using two additional datasets of thermal and multispectral images, to demonstrate our multi-sensor spatial co-registration methodology. Input data were collected at three study sites in the Windmill Islands region, Antarctica (near the Australian base, Casey), where some of the most well-developed continental Antarctic vegetation is located.

The three study sites were named: Robinson Ridge, Red Shed, and Antarctic Specially Protected Area 135 (ASPA 135) (see Figure 1 for location overview). Robinson Ridge is located approximately 10 km south of Casey station, the Red Shed site is beside a melt lake behind the main accommodation building at Casey, and ASPA135 lies approximately 500 m east of Casey. Further description of all three test sites can be found in Turner et al. [26] and in Lucieer et al. [25].

2.2. Platform

Multi-rotor UAVs are becoming more commonplace and are frequently used for commercial and recreational aerial photography. For this study we used an eight rotor Mikrokopter Micro-UAV called an “Oktokopter” supplied by HiSystems (GmbH, www.mikrokopter.com, Germany) (see Figure 2). The Oktokopter had a payload capacity of around 1 kg, a flight duration of 5 min (with a typical payload), and was equipped with a gimballed camera mount (i.e., self leveled during flight based on onboard gyroscopes that measure the roll and pitch of the airframe), to which we individually fitted each of the three sensors. The Mikrokopter flight electronic systems can be used to automatically maintain level flight, control the altitude, log system data, and to fly the UAV through a series of predefined, three-dimensional waypoints.

2.3. Visible Digital Camera

To collect visible imagery we used a Canon 550D Digital Single Lens Reflex (DSLR) camera (15 Megapixel, 5184 × 3456 pixels, with Canon EF-S 18-55 mm F/3.5-5.6 IS lens). The image capture rate was controlled by the UAV’s flight control board, which was programmed to emit a pulse at a desired frequency. The flight control board was connected to a custom-made cable that triggered the remote shutter release of the camera. The Canon camera was operated in shutter priority mode (a fast shutter speed was required to minimize motion blur), in which the desired shutter speed (typically 1/1250–1/1600 sec) was set before flight and the exposure was adjusted automatically by varying the camera’s aperture. Images were captured in RAW format and stored on the memory card in the camera for post-flight download.

2.4. Thermal Infrared Sensor

To collect TIR imagery we used a FLIR Photon 320 (FLIR Systems, Inc. the USA, www.flir.com) uncooled thermal sensor (see Figure 2). The Photon 320 had a 14 mm lens providing a 46° field of view and acquired image frames of 324 × 256 pixels as raw 14-bit Digital Numbers (DNs) at the rate of 9 Hz. Image frames from the camera were converted into ethernet data packets by the FLIR Ethernet module and this data was then stored on a Single Board Computer (SBC), a Gumstix Verdex Pro XM4-BT, equipped with netCF and console expansion cards. System time of the Gumstix SBC was set to GPS time prior to flight, so that the thermal data files could be synchronized with UAV GPS log files. The 9 Hz data rate was too fast for the data buses of the Gumstix SBC and thus images could only be collected at a rate of around 1 Hz, which was fast enough for our purposes and retained sufficient image overlap on ground.

After the flight, the raw image data was downloaded from the SBC memory card and processed with code written in the IDL/ENVI image-processing environment (Exelis Visual Information Solutions, Inc. USA, www.exelisvis.com) to extract the image frames from the captured data packets. Extracted images were stored as 16-bit ENVI single band files containing the original 14-bit raw DNs as collected by the TIR sensor. A set of JPEG quick look images were simultaneously generated allowing a visual check of image quality. As the DNs typically do not cover the full 14-bit dynamic range, a contrast stretch to the data had to be applied, so that subsequent image processing software was able to identify features within the images. To identify an appropriate stretch, we created a histogram of the DNs of all pixels in all images, chose an upper and lower threshold such that the full dynamic range of the scene was covered, and applied a linear stretch based on these minimum and maximum thresholds. The images were then stored as 16-bit TIFF files, keeping note of the thresholds used such that the pixel values could later be converted back to the original DNs. The DN-values in the thermal imagery represent at-sensor radiance. After mosaicking and co-registration the DN-values were converted to absolute temperature in °C based on an empirical line correction. Nineteen targets with similar emissivity (0.97 assumed for moss and dark rock) were marked with shiny aluminum disks. Due to the very low emissivity of shiny aluminum these targets were clearly visible in the thermal imagery. A temperature observation was collected with a thermal radiance gun (Digitech QM7226) in between two aluminum disks. A GPS coordinate was also recorded for these observations. Matching pixels were extracted from the thermal imagery and based on the matching reference temperatures a linear regression was calculated. With this empirical relationship we converted the whole thermal mosaic into absolute temperature, assuming a constant emissivity of 0.97 (which is justifiable given our interest in the moss bed).

2.5. Multispectral Sensor

The multispectral sensor used in our study was a Tetracam (Tetracam, Inc. the USA) mini-MCA (Multiple Camera Array) with an array of six individual image channels. Each channel has its own Complementary Metal Oxide Semiconductor (CMOS) sensor that could acquire 10-bit image data at an image size of 1280 × 1024 pixels. It was possible to fit customized waveband pass filters to each lens, allowing the user to define the desirable spectral band configuration. The data sets collected for this study had 530, 550, 570, 670, 700 and 800 nm optical filters fitted, with a Full-Width at Half Maximum (FWHM) of 10 nm.

The mini-MCA could be set to a “burst” mode in which it captured images continually from the time the shutter release was sfirst pressed, we used the maximum rate of 0.5 Hz. Each of the six arrays stored the images in a proprietary raw format onto individual compact flash memory cards. After the flight the image data was downloaded from the cards resulting in six files for each camera exposure. As with the TIR data, we designed processing code in the IDL/ENVI environment to read the raw format files and to merge the layers into a single six band, 16-bit (to store the 10-bit data) ENVI image file format. The next stage was to correct mini-MCA imagery for sensor noise and other image distortions. Detailed description of these corrections can be found in Kelcey and Lucieer [27]. In short, the following three corrections were applied: (i) noise reduction using dark current imagery, (ii) lens vignetting correction based on spatially dependent correction factors, and (iii) a Brown-Conrady model removed lens distortion [27].

Finally, the image bands had to be aligned, as the six mini-MCA camera lenses were spatially offset. Tetracam Inc. provides software and alignment equations to correct for these offsets, but similar to Laliberte et al. [11] we found this alignment correction inaccurate at our typical UAV flying heights. To improve the band misalignment, we developed our own technique based on detecting geometric features within the imagery with the Scale Invariant Feature Transform (SIFT) keypoint detector [28]. The alignment process considered the first spectral band (Band 1) of each image to be the master and aligned the other bands to it by matching key points between the bands. SIFT was run on each band of an image to create a set of key files and a key matching algorithm was then run between Band 1 and each of the other bands. Extracting the x, y locations of the matching features allowed us to create a control point file that aligned a given band with Band 1, i.e., there were five control point files for bands 2 to 6. The control points were applied to each band with a Delaunay triangulation combined with a nearest neighbor resampling in order to perform a non-linear local transformation for alignment with the first band. The band alignment was highly dependent on the distance between the camera and the imaged surface. Provided that this distance remained relatively constant (within ±5 m), the technique would create a set of alignment parameters applicable to the current dataset. Final inspection of the band alignment for objects with sharp edges in the images was found to be satisfactory. Moreover, the method had the added advantage of being fully automated.

2.6. Mosaicking of Visible Imagery

During a typical flight of our UAV the Canon camera collected around 200 images. It was therefore necessary to join the images into a single mosaic of the whole study area. The image mosaic also needed to be georeferenced, such that the imagery from the different sensors could be co-registered. There are various methods for mosaicking UAV imagery; e.g., Berni et al. [8], Bryson et al. [29], Laliberte et al. [30], Turner et al. [26], and Turner et al. [31]. Recently, new commercial software packages for automatically georectifing UAV imagery have become available (a review of some of these packages can be found in Turner et al. [31]). Based on our previous research results we selected Photoscan Professional by Agisoft (Agisoft LLC, Russia) to georectify and mosaic the visible UAV imagery.

An overview of the Structure-from-Motion (SfM) workflow in Photoscan software can be found in Lucieer et al. [25]. Prior to processing, the images were geotagged with their approximate location as recorded by the UAV’s on-board navigation-grade GPS. The internal time of the camera was set to GPS time prior to flight to ensure that the images could be easily synchronized with the position data in the UAV GPS log file. The blurry images were then detected with an algorithm that calculates an image blur metric according to Crete et al. [32]. A more detailed description of this image processing stage can be found in Turner et al. [31]. In the final step images with excessive overlap are removed [31].

Once the set of images to be processed was finalized, they were imported into Photoscan, which then detected and matched thousands of features between the images. Using these matches it performed a bundle adjustment to estimate the camera positions, orientations, and lens calibration parameters. Based on this information, the geometry of the scene (in the form of a 3D model) was created by applying a dense, multiview stereo reconstruction to the aligned images. Once the 3D geometry of the scene was constructed, a Digital Surface Model (DSM) and an orthophoto mosaic could be exported [25].

To improve the absolute spatial accuracy of the mosaics, we manually located Ground Control Points (GCPs) distributed within the imagery. The GCPs were 30 cm diameter metal disks with a bright orange rim that were laid out in the study area prior to flight and measured with Differential RTK GPS (DGPS) with a typical accuracy of 2 cm in the horizontal and 4 cm in the vertical direction (relative to a local coordinated benchmark). Similar to the large GCPs, between 20 and 45 smaller (10 cm diameter) orange metal disks were randomly laid out across the study area and also coordinated with a DGPS. These small GCPs were used later as check points to verify the accuracy of the georeferenced mosaics.

Photoscan provided a simple interface to mark the location of a GCP on the 3D model and its location was then automatically marked on all the images that covered that part of the model. The user then needed to manually verify and, if necessary, adjust the location of the GCP in each image. Although this process was time-consuming, taking from 2 to 4 h to mark about 20 GCPs in a dataset of 200 photographs, the significant improvement in spatial accuracy justified the work required.

2.7. Mosaicking of Thermal Infrared Images

Similarly to the visible camera, the TIR sensor could collect hundreds of images in a single flight and it also required selection of the best images from the acquired dataset. The response time of the microbolometer in the FLIR Photon 320 is approximately 10 milliseconds (see www.flir.com), giving it an effective shutter speed of 1/100th of a second, resulting in motion blur in around 40% of the images. Blurriness was assessed with the algorithm described in Section 2.6 and images with a blur metric greater than 0.3 were automatically removed from the dataset (see Figure 3). TIR images were subsequently processed with Photoscan in a manner similar to the visible imagery.

Photoscan was provided with initial estimations of camera position based on the time stamp of the image and the position as recorded in the on-board GPS log file. The TIR imagery was stored as a single band 16-bit file. Features in the images were enhanced with a linear contrast stretch, which was based on the minimum and maximum temperature DNs detected in the full scene (see Section 2.4). Once Photoscan had aligned the images, the GCPs needed to be identified within the imagery to allow accurate co-registration. It should be noted that only large aluminum trays (30 cm diameter) with an unpainted central part (see Figure 4b) were detectable in the lower resolution TIR imagery (10 cm/pixel for a typical UAV flight). The metallic aluminum surface of the trays with a low emissivity appeared as very cold (dark) pixels, which made the GCPs easy to visually identify. After the GCPs were marked, the scene’s geometry was constructed, and an orthophoto was generated following the same approach as described in Section 2.6.

2.8. Mosaicking of Multispectral Images

Unlike the TIR and the visible imagery it was not possible to process the multispectral imagery with Photoscan. Although multi-band 16-bit images could be imported into Photoscan, it was not feasible to reliably align the images of different spectral bands, no matter what parameters were selected. The CMOS sensors in the mini-MCA have a rolling shutter, which built up each image as a scan from top to bottom rather than a whole-frame snapshot as in the case of a global shutter. Given the movements of the sensor during image acquisition the rolling shutter leads to geometric distortions in each image, which were unpredictable. The SfM algorithm in Photoscan expects images to be acquired by a global shutter. Thus, the distorted image geometry of the mini-MCA leads to very poor and false image matching results. Laliberte et al. [33] developed a method for rectifying, georeferencing, and mosaicking UAV visible imagery by matching the individual images with a pre-existing orthophoto of the study area using image correlation techniques. This method served as the basis for our approach, however, we have implemented the SIFT algorithm as we did with the mini-MCA band alignment workflow (see Section 2.5). Instead of using a low-resolution orthophoto we matched each mini-MCA frame to the ultra-high spatial resolution visible orthophotomap (Section 2.6). This allowed SIFT to select thousands of features per image to be used as control points, providing a denser transformation matrix.

The mini-MCA collected a large number of overlapping images during flight. High overlap (80%–90%) is essential for the SfM algorithm, but for the mini-MCA imagery we only needed overlap of about 25%–30% to form a continuous mosaic. Therefore, the first processing step of the multispectral mosaicking workflow was determining an optimal subset of input images. The position of the airframe at the time of each image exposure, logged by the on-board GPS unit, is used to select the mini-MCA images based on their spatial distribution. Through process optimization we determined that a threshold of 7 m between the image positions achieved a mosaic with full coverage, but minimal seam lines. This selection reduced the number of images in a test dataset from 142 to 41, greatly reducing the number of seam lines within the final mosaic, whilst maintaining sufficient image overlap of around 25% (as determined via experimentation across all datasets). The blurriness of mini-MCA images was verified in the same manner as the other two datasets. However, the blur factor of the mini-MCA images was, in general, so low that no images were excluded from the datasets.

Before we matched features in a mini-MCA image, we had to identify all the features in the visible orthophoto. To overcome limitations of the conventional SIFT algorithm implementation that usually runs only on low-resolution images, we implemented a SIFT distribution called “libsiftfast” [34], which uses multiple CPU cores and can process large images. As result, SIFT was able to detect 3,643,780 features in one of our 8690 × 17,215 pixels large visible mosaics. Since the SIFT algorithm runs per single image band, the six mini-MCA spectral bands supply more geometrical features for the matching process than the standard RGB imagery. From the millions of matches, there are often many false matches, which were removed with a Random Sample Consensus (RANSAC) algorithm developed by Fischler and Bolles [35]. In the application of the RANSAC algorithm we made the assumption that we were working with a projective model and that the epipolar constraints would hold despite the possibility of rolling shutter distortions within the imagery. To ensure these distortions did not affect the results we applied a low distance threshold (0.01 in normalized space), which reduced the number of resultant matches, whilst ensuring there were no false positives. We did not require a large number of matches as the remaining feature matches (post RANSAC) were not going to be used in an SfM context but were instead used to create control point files.

Using the feature matches, and the feature key files containing the locations of each feature in each image, we automatically created a series of dense GCP files for each mini-MCA image. This was facilitated by matching the known x, y coordinates of mini-MCA image features with their corresponding x, y coordinates in the RGB orthophoto, which were converted to easting and northing coordinates as the RGB mosaic was already georeferenced. The GCPs were then used in a Delaunay triangulation, which transformed a mini-MCA image using a nearest neighbor resampling algorithm into the same coordinate system as the visible and TIR mosaics. The pixel size of the resulting georeferenced image was set based on the flying height during the image acquisition (3 cm/pixel for our test datasets). Finally, all the images were merged into a georeferenced mosaic with the IDL/ENVI mosaicking routine using a feathering of the seam lines to smooth the original image edges.

2.9. Calculating MTVI2 from Multispectral Data

Raw images of the mini-MCA, recorded as 10-bit digital counts with maximal signal strength equal to 1024 DNs, were transformed into the physically meaningful relative reflectance by applying an empirical line correction [36]. Three spectrally flat calibration panels (40 × 40 cm) of white, grey, and black color with reflectance intensities ranging from 2 up to 80% were placed within the UAV flight path and captured in the mini-MCA imagery. The target’s actual DN values (approx. 130–140 pixels per panel) were extracted from the airborne image and empirically related to their reflectance functions measured with a spectrally calibrated ASD HandHeld2 (HH2) spectroradiometer (Analytical Spectral Devices, PANalytical Boulder USA) on the ground immediately after completion of the mini-MCA acquisition. The absolute reflectance of the calibration targets measured with the ASD-HH2 between 325 and 1075 nm (751 spectral bands of 1 nm bandwidth) were spectrally convolved to resemble the six broader spectral bands of the mini-MCA instrument with FWHM of 10 nm. The empirical line correction coefficients established between the convoluted ASD-HH2 and the six acquired mini-MCA spectral bands were then applied per pixel to the mini-MCA image mosaic to remove optical attenuation caused by scattering and absorption processes of atmospheric gases and aerosols between the sensor and observed surfaces and to standardize the multispectral signal as the relative reflectance function.

An efficient way of detecting photosynthetically active vegetation in multispectral imagery and assessing its actual physiological state is to transform the relative reflectance function into an optical vegetation index. The Modified Triangular Vegetation Index 2 (MTVI2) was originally introduced by Haboudane et al. [37] to estimate green biomass density of spatially homogeneous agricultural crops. It is computed as:

Being stressed by insufficient water supply and high photosynthetically active and ultraviolet irradiation, Antarctic moss turf can, within few days, change in compactness and pigmentation. Over days to weeks it changes from a healthy, green open-leaved form to a stress-resisting denser, yellow-brown or red closed packed turf. If dry periods persist over longer periods (months to years) the turf will lose photosynthetic pigments, forming grey-black mounds of moribund (dormant/dead) moss. Since these physiological stress reactions systematically alter the moss reflectance at the wavelengths of 550, 670 and 800 nm, we could apply MTVI2 to separate photosynthetically active moss (health > 60%) from moribund moss, lichens, and the rocky surroundings.

For assessment of the ability of the MTVI2 index to determine moss health, the Robinson Ridge MTVI2 results were compared with field samples from the 2012 field season collected as part of a long term monitoring system for Australian State of the Environment Indicator 72 (SoE 72). Established in 2003, this monitoring system comprises a set of 30 permanent quadrat locations, across 10 transects spanning a water gradient across three community types: from the wettest community dominated by mosses (Bryophyte community) to the driest community dominated by moribund moss encrusted with lichens (Lichen community), with a Transitional community between [39]. For each quadrat, the percentage of live bryophytes was evaluated from 9 small samples, each containing approximately 20–50 moss shoots, at 9 intersections within a 20 × 20 cm grid. For a more detailed description of moss ground monitoring see Wasley et al. [39] and Lucieer et al. [25] and for a description of spectral properites of Antarctic moss see Lovelock and Robinson [40]. To replicate the ground based sampling scheme, a grid of 3 × 3 pixels was extracted from the MTVI2 map and averaged per sampled quadrat.

2.10. Accuracy Assessment

To achieve the best possible co-registration of the three different datasets, it was essential that each of the mosaics was georeferenced with the highest spatial accuracy possible. For the RGB orthomosaics we, therefore, measured the positional error of all small orange disk check points that were measured with DGPS, but not used by Photoscan to georeference the mosaics. These disks were, unfortunately, too small to be visible in the TIR imagery and to be accurately identified in the mini-MCA imagery. For these two datasets the larger GCPs were used as check points. The Root Mean Square Error (RMSE) was computed between check point coordinates measured in the field with DGPS and coordinates retrieved from georeferenced image mosaics to assess the overall spatial accuracy of each dataset.

3.1. Image Mosaics

The three tests sites were flown with our manually navigated Micro-UAV on three separate days. Manual navigation was required due to failure of autopilot navigation, which was caused by the extreme magnetic declination of Eastern Antarctica (∼100° West). A basic description of the acquired datasets is provided in Table 1. Flights were all carried out in good solar illumination conditions, light winds, and at an altitude of approximately 50 m Above Ground Level (AGL). Despite flying at the same height AGL during each UAV mission, the resulting image mosaics have different spatial extents and different spatial resolutions caused by differences in the technical parameters of each sensor.

Applying the methods described in Section 2, we created georeferenced mosaics for each of the nine flights. An example of a visible mosaic from the Robinson Ridge site is shown in Figure 4a. Figure 4b illustrates the detail that can be seen in this RGB imagery with a pixel size of 1 cm. The same spatial subsets of the lower resolution thermal infrared mosaic (10 cm/pixel) and of the false colour multispectral mosaic (3 cm/pixel) are shown in Figure 4c,d, and finally Figure 4e gives an example of a typical spectral signature for the healthy moss.

Using the methods described in Section 2.10, the spatial accuracy of each orthomosaic was measured and a summary of the spatial errors and RMSE is provided in Table 2. The RMSE ranges in general from 1 to 2.6 pixels, which means that all the mosaics exhibit a comparable level of spatial accuracy. The ASPA135 visible mosaic has the lowest RMSE, which can be attributed to the fact that it is a relatively small site with limited geomorphological variability, whereas the other two test sites are larger in spatial extent and have more diverse terrain morphology.

3.2. Co-Registration Accuracy

Large GCPs were used as cross-comparison check points for the co-registration accuracy assessment, since they were the only features clearly identifiable in all three datasets. The location of the comparison points in the visible mosaic was considered to be the reference position. The co-registration errors of these GCPs in the mini-MCA and TIR mosaics are summarized in Table 3. The RMSEs are generally around 2 pixels, which matches the absolute RMSE that was obtained for all the mosaics listed in Table 2.

3.3. Assessing Health of Antarctic Moss from Multisprectral Imagery

To demonstrate sufficient spatial accuracy of these geocoded mosaics and the ability of the multispectral imagery to assess the actual health state of Antarctic mosses, we computed the MTVI2 optical vegetation index from the mini-MCA imagery collected at the Robinson Ridge test site (see Section 2.9). The Robinson Ridge data was selected as the most suitable of the three available datasets, because it has the least snow cover, it was acquired under optimal light conditions (bright, but diffuse irradiation), and ground observations of actual moss health are available for this site. Since reflectance signatures of stressed mosses and agro-systems with low leaf density (i.e., low leaf area index) are spectrally similar, we could apply MTVI2 to assess the spatial distribution of the health state of moss bed at Robinson Ridge.

Comparing the moss health measured in summer 2012 and mean MTVI2 values derived from the mini-MCA mosaic we found a strong quadratic relationship between the two (see Figure 5). MTVI2 is insensitive to moss health in the driest lichen community quadrats due to the low abundance of photosynthetically active moss. MTVI2 does, however, show a statistically significant positive relationship for quadrats with a significant and/or dominant presence of bryophytes (R² = 0.636). Such a strong correlation provides evidence of a good co-registration agreement between mini-MCA and DGPS localization of quadrats and georectification of the mini-MCA image mosaic.

The statistical relationship shown in Figure 5 allowed us to approximate per-pixel moss health of the whole moss bed captured in the mini-MCA mosaic. Unfortunately, a considerable number of erroneous pixels were identified after a close inspection of the moss health map. These were caused, for instance, by high sensor noise combined with low light in shadows that produced incorrect spectral signatures of rock surfaces, mimicking high MTVI2 values of healthy moss. To prevent false moss health estimates, we applied the following rules: (i) the calculated per-pixel value of the Normalized Difference Vegetation Index (NDVI) must always be positive, (ii) the reflectance value of the 800 nm spectral band must be positive and greater than the reflectance in the shorter wavelengths, i.e., at 550, 570 and 700 nm, and finally (iii) to eliminate negative estimates, moss health is assessed only if the MTVI2 value is greater than 0.4. The final map of moss health with these rules applied is displayed in Figure 6b.

3.4. Assessing Temperature of Healthy Moss from Thermal Infrared Image Msaic

Simultaneously with the UAV flight over the Robinson Ridge study site, ground measurements of surface temperature were carried out with a handheld infrared thermometer at various locations recorded with DGPS. The corresponding DN value at each of the ground sample points was extracted from the georeferenced TIR mosaic and compared to the ground temperature measurements. A strong linear relationship established between both datasets was used to convert the TIR DN values to indicative surface temperatures (see Figure 7), providing a map of indicative surface temperature for the whole study site. Applying the moss health map as a mask for the TIR imagery produced the surface temperature map of only moss pixels. The resulting surface temperature map of moss with a relative health score greater than 60% is shown in Figure 6c.