CORRECTION OF AIRBORNE PUSHBROOM IMAGES ORIENTATION USING BUNDLE ADJUSTMENT OF FRAME IMAGES

To compute hyperspectral orthophotos of an area, one may proceed like for standard RGB orthophotos : equip an aircraft or a drone with the appropriate camera, a GPS and an Inertial Measurement Unit (IMU). The position and attitude data from the navigation sensors, together with the collected images, can be input to a bundle adjustment which refines the estimation of the parameters and allows to create 3D models or orthophotos of the scene. But most of the hyperspectral cameras are pushbrooms sensors : they acquire lines of pixels. The bundle adjustment identifies tie points (using their 2D neighbourhoods) between different images to stitch them together. This is impossible when the input images are lines. To get around this problem, we propose a method that can be used when both a frame RGB camera and a hyperspectral pushbroom camera are used during the same flight. We first use the bundle adjustment theory to obtain corrected navigation parameters for the RGB camera. Then, assuming a small boresight between the RGB camera and the navigation sensors, we can estimate this boresight as well as the corrected position and attitude parameters for the navigation sensors. Finally, supposing that the boresight between these sensors and the pushbroom camera is constant during the flight, we can retrieve it by matching manually corresponding pairs of points between the current projection and a reference. Comparison between the direct georeferencing and the georeferencing with our method on three flights performed during the Leman-Baikal project shows great improvement of the ground accuracy. 1. CURRENT METHODS FOR PUSHBROOM DATA REGISTRATION 1.1 Pushbroom Sensors Pushbroom sensors are common to acquire hyperspectral data as they usually propose high spatial resolution and a high number of bands. They work as line sensors, and they are usually used on moving platforms to scan zones (see Figure 1).


Pushbroom Sensors
Pushbroom sensors are common to acquire hyperspectral data as they usually propose high spatial resolution and a high number of bands. They work as line sensors, and they are usually used on moving platforms to scan zones (see Figure 1).
optical center x y z Figure 1. Pushbroom sensors are optical sensors consisting in a line of pixels. They are usually moved orthogonally to the line to construct images.
To establish a hyperspectral cartography of an area of interest, it is possible to mount such a sensor on an aerial vehicle and scan the area by flying over it. Each of the acquired hyperspectral lines has six orientation parameters, called Exterior Orientation Parameters (EOP): three coordinates representing its position, and three angles representing its attitude. These parameters are relative to a coordinate system; in the following, we work with the WGS84 Geoid. The position parameters we use are latitude, longitude and altitude, and the attitude parameters are roll, pitch and yaw angles, defined as in Figure 2. Inherent parameters of the camera, called Interior Orientation Parameters (IOP) also impact the georeferencing task: the focal length of the camera, its principal point location and the radial distortion of its lens. When all the orientation parameters are known for a single line, this line can be georeferenced.

Georeferencing
One can achieve decent georeferencing of airborne imagery by getting these parameters through integration and Kalman filtering of data provided by a GPS and an IMU (see, for instance, (Mostafa and Hutton, 2001)). This process is called Direct Georeferencing. However, errors remain, due to the accuracy of the navigation sensors, the unperfect synchronization between them and the cameras, or IOP not known precisely enough. It is especially critical for pushbroom sensors where a coherent ortho-rectified image can be obtained only if the EOP of each of the acquired hyperspectral lines are known accurately. Figure  3 shows a projection of hyperspectral data where the navigation data used is obtained by integrating the information provided by the GPS+IMU system Ekinox-N designed by SBG Systems. For a single flight line, trajectory models using random processes allow to achieve better ground precision (see (Lee and Bethel, 2001)), but in the general case, post-processing, using Ground Control Points (GCP) or extra information, is needed to compute a coherent and well georeferenced image. This is called Indirect Georeferencing. The GCP approach consists in matching points in our data with coordinates in the reference (i.e, on the Earth when these points have been measured on the field) and then applying corrections to the navigation data according to this matching. In the case of a pushbroom sensor, one would need 3 GCPs per line to determine the 6 orientation parameters of a given line (see (Lee et al., 2000)), which is not feasible. But combining the direct and the indirect approach allows good georeferencing: a few GCPs can help to geometrically correct data after a first raw projection, using interpolation methods like Particle Swarm Optimization (see (Reguera-Salgado and Martín-Herrero, 2012)).

OUR APPROACH
In this section, our approach is detailed. It is decomposed as follows: we perform a bundle adjustment of the RGB images captured by our RGB camera during the flight. It utilizes the navigation data output by integrating the GPS/IMU information. Then, the corrected navigation data for the RGB camera is retrieved, and used to estimate the boresight between this camera and the IMU, as well as the corrected navigation data for the IMU, using a leastsquare compensation. The last step consists in projecting the hyperspectral data using this new navigation data and performing a last correction, taking into account the boresight between the pushbroom sensor and the IMU; this boresight is estimated using tie points between the latest projection and the reference. Figure  4 summarizes this process.

Bundle Adjustment of RGB Images
Using feature matching techniques, RGB images can be stitched together to create a 3D model of the area. If a priori navigation data is also given, the process can correct it as well as the IOP like the focal length of the camera, the position of its principal point and the radial distortion. This process is called bundle adjustment (see (Triggs et al., 2000) for reference). In our case, we use the navigation data integrated by the Kalman filter. Even though it is different from the navigation data for the RGB camera (since there is a boresight and a lever-arm between the two), it is still close enough to be considered as a good first estimation of the orientation parameters of the camera. The bundle adjustment is processed using Agisoft PhotoScan (see (Agisoft, 2013)); Figure 5 shows a 3D model ofÉcole Polytechnique Fédérale de Lausanne (EPFL) produced by this software. From there, we can retrieve corrected navigation data for the RGB camera.

Navigation Sensors Data Correction
The orientation parameters output by the previous step are for the frame camera. The correct parameters for navigation sensors are not the same, since there is a lever-arm as well as a boresight between them and the camera. Since our hyperspectral camera has a frequency 50 times higher than the one of our RGB camera, only one pushbroom line out of 50 can be directly corrected; we need to use the data acquired by navigation sensors to interpolate between the corrections.

Position Correction
The position of the IMU will not be computed. Indeed, in our system, the norm of the lever-arm between the pushbroom camera and the frame camera is less than 10 cm. Therefore, we use the position parameters retrieved from the bundle adjustment step for the pushbroom camera. As The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B3, 2016 XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic a consequence, navigation sensors data can be used as such to interpolate between the corrections of the camera position.
Consider a position parameter p among latitude, longitude and altitude. Let pr and pc be respectively the raw value of this parameter (output of the Kalman Filter) and the corrected value. Let t0 be the timestamp of a line corrected thanks to the bundle adjustment. So is the line which timestamp is t0 + 50∆t, where ∆t is the sampling time of the hyperspectral camera. Since the IMU is a drifting sensor (accurate for short periods of time, but getting less and less accurate as time goes by), the variations it gives around t0 are highly reliable close to t0, and less when getting closer to t0 + 50∆t. Same goes for the variations around t0 + 50∆t. We therefore interpolate pc, for t ∈ [t0, t0 + 50∆t], using Equation 1.
is the weight given to the navigation data close to t0. Its value is given by Equation 2. (2)

Attitude Correction
The attitude of the navigation sensors has to be computed. Indeed, to interpolate the attitude data, it is not possible to use the data from the navigation sensors as such, since it refers to its own frame, while the attitude given by the bundle adjustement refers to the RGB camera frame. At a given location, let R Camera N ED be the rotation matrix from the North-East-Down local-level frame to the camera frame, R IM U N ED the one from the local-level frame to the IMU and R bore = R Camera IM U the rotation representing the boresight from the IMU to the camera. The relation between these matrices is given by Equation 3.
Each of these rotation operators can be characterised by its roll, its pitch and its yaw, applied in this order: roll rotation around the x-axis, pitch rotation around the y-axis, and yaw rotation around the z-axis. Let r, p, y be these 3 parameters; the corresponding rotation operator is given by Equation 4. R(r, p, y) =   cos y cos p cos y sin p sin r−sin y cos r cos y sin p cos r+sin y sin r sin y cos p sin y sin p sin r+cos y cos r sin y sin p cos r−cos y sin r − sin p cos p sin r cos p cos r   This boresight is constant during the flight as the sensors are rigidly attached to the vehicle. We call (r b , p b , y b ) its roll, pitch and yaw. Similarly, at time t, let (rn(t), pn(t), yn(t)) be the orientation parameters output by the navigation sensors, and (rc(t), pc(t), yc(t)) the one for the RGB camera (given by the bundle adjustment). Let k be the number of RGB images available, and (ti), i ∈ [1, k] the times of acquisition of each image.
The drifts (errors) in the orientation parameters added to the current drifts at time ti are called (dri, dpi, dyi). For each time ti, i ∈ [1, k], we can rewrite Equation 3 using the orientation parameters of each rotation matrix. The set of equations obtained is given by 5.
At time ti, (dri, dpi, dyi) are small, unlike the total drifts ( i j=1 drj, i j=1 dpj, i j=1 dyj), where all the drifting errors and the subsequent noise from gyrometers have been added. Therefore, this model allows to keep the time-dependency of the drifts while having parameters that shall be small. Boresight parameters shall also have low values, as the sensors are rigidly attached, pointing in the same direction, in our system. By retrieving the orientation parameters from the left-hand side and the right-hand side of each equation in 5, we get a set of 3k equations with 3k + 3 unknown parameters (all the drifts plus the 3 parameters of the boresight). This can be seen as an optimization problem: consider the vector of all the parameters to be determined, as given by Equation 6.
Our goal is to minimize v T v under the constraints given by the set of equations 5. With all parameters equal to 0, the equations are not verified, and for each equation there is a gap between the lefthand side and the right-hand side. If we call w the column vector which elements are all these gaps and B the Jacobian matrix of the set, then it is known, using Lagrange multipliers theory, that v can be approximated by Equation 7.
As a result, we get all corrected parameters for the navigation sensors for the times of RGB acquisitions, as well as a side-product The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B3, 2016 XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic (the boresight). The rest of the corrected navigation sensors data is interpolated using the exact same interpolation method as described in 2.2.1.

Pushbroom Camera Attitude Correction
Once the attitude from the navigation sensors is known, the last step is to find the boresight between them and the pushbroom camera. A georeferenced RGB orthophoto and a georeferenced digital elevation model can be imported from Agisoft PhotoScan to our software HypOS, based on the NASA Worldwind SDK; overlooking at this unknown boresight, a first projection of the hyperspectral data can be made. Then, we manually point at corresponding pairs of points between the RGB data and the hyperspectral data. This is possible as the hyperspectral data can be seen as RGB, using 3 bands of the spectrum roughly corresponding to red, green and blue. With a large enough set of pairs, it is then possible to find the best boresight, as well as some IOP (focal length, principal point coordinates), to match these pairs. In the Earth-centered Earth-fixed frame, consider P P = (xpp, ypp, zpp) as the principal point of the pushbroom camera, and G = (xG, yG, zG) as a ground point we want to match; let (r, p, y) be the three orientation parameters of the boresight from the IMU to the camera. In the front-right-down frame of the camera centered at the principal point, the coordinates of a pixel are (δx, ypix + δy, f + δf ) where ypix is the coordinate of the pixel along the line of the pushbroom and f is our input estimation of the focal length. If the principal point was centered relatively to the camera, δx and δy would be zero; but it is not rigourously true, therefore δx and δy are interior orientation parameters that we might want to determine accurately. For a given pair, Equation 8, also known as the collinearity equation, should be verified.
F (r, p, y, δx, δy, δf ) = G − P P G − P P − R(r, p, y)(δx, δy + ypix, f + δf ) R(r, p, y)(δx, δy + ypix, f + δf ) = 0 (8) 6 unknown parameters have to be determined. This is possible as long as 6 pairs of points are input. Then, these parameters can be determined by a least-squares approximation. Let k be the number of input pairs, (Fi), i ∈ [1, k] the k corresponding functions as described in Equation 8, v the vector of the six parameters to be found, and w the vector of the gaps (Fi(0, 0, 0, 0, 0, 0)), i ∈ [1, k]. The estimation of the parameters can be refined using Equation 9.
This compensation can be applied iteratively till the correction is smaller than a given threshold. The output includes the boresight and the correction of interior orientation parameters; with these, navigation data for the pushbroom sensor is computed, and then we can georeference the data.

RESULTS
To assess the performance of our method, we compare here the precision of the projection before and after using it, on 3 flights of different mean elevations. These flights were performed in the context of the Leman-Baikal project, which consists in a study of both lake Geneva (Switzerland) and lake Baikal (Russia) by making aerial surveys with hyperspectral cameras (see (Akhtman et al., 2014)). The cameras, together with the navigation sensors, are mounted on a ultralight trike. We have observed an important yaw drift along our flights, as well as a boresight, which resulted in bad georeferencing of the data. Therefore, this data is a good case study for our algorithm.
Our reference can either be the orthorectified photo output by PhotoScan or Bing Maps Imagery. The quality of each projection is characterized by its mean residual over a set of 50 reference ground control points. This has been tested on flights of different mean elevations. Results are summarized in Table 1 Quantitatively, images produced using our algorithm are much better georeferenced than the one created using Direct Georeferencing. The correction gets more important with the altitude: for flights with a higher elevation, the mean residual using Direct Georeferencing grows up to 50 meters, while the impact of the height is way less important with our method. The Ground Sampling Distance (GSD) of our system is roughly a thousandth of the height. Consequently, the order of the residual is 1-6 GSD for the lower altitude flight, and 1-4 GSD for the other two. Qualitatively, Figures 6, 7 and 8 show, for each of these flights, the flight path, the raw projection (before correction) and the projection after correction. The projected data is displayed with a lighter color to differentiate it from the reference (from Bing Maps Imagery) in the background. Particular features -roads, buildings or rivers -are highlighted in specific colors. They show that some shifts or distortions (due to misprojected overlapping lines) have been corrected in the processed images, which fit much better the background reference images and are more coherent.

CONCLUSION
We have presented a method for correcting geometrically hyperspectral pushbroom data: the raw projection is rectified using the data retrieved from the bundle adjustment of RGB images. Then, least-squares approximation allows to determine precisely the boresights between the sensors embedded in the vehicle and to retrieve the orientation parameters of each hyperspectral line to perform coherent projection on a reference globe. Quantitatively, this process drastically reduces the average distance error between projected points and their corresponding ground control points, compared with the direct georeferencing. This method can be used to georeference any airborne data acquired together with navigation data and RGB data, even with large attitude variations within the flight. In our future work, we will be interested in automating the matching between the reference image and the data to be ortho-rectified. Using image processing techniques, we can hope to perform feature matching between the two sets of data and make the algorithm fully unsupervised.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B3, 2016 XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B3, 2016 XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic