Bundle Adjustment Using Space-Based Triangulation Method for Improving the Landsat Global Ground Reference

There is an ever-increasing interest and need for accurate georegistration of remotely sensed data products to a common global geometric reference. Although georegistration has improved substantially in the last decade, the lack of an accurate global ground reference dataset poses serious issues for data providers seeking to make geometrically stackable analysis-ready data. The existing Global Land Survey 2000 (GLS2000) dataset derived from Landsat 7 images provides global coverage and can be used as a reference dataset, but its accuracy is much lower than what can be attained using the agile and precise pointing capability of the new spacecrafts. The improved position and pointing knowledge of the new spacecrafts such as Landsat 8 can be used to improve the accuracy of the existing global ground control points using a space-based triangulation method. This paper discusses the theoretical basis, formulation, and application of the space-based triangulation method at a continental scale to improve the accuracy of the GLS-derived ground control points. Our triangulation method involves adjusting the spacecraft position, velocity, attitude, attitude rate, and ground control point locations, iteratively, by linearizing the non-linear viewing geometry, such that the residual errors in the measured image points are minimized. The complexity of the numerical inversion and processing is dealt with in our approach by processing and eliminating the ground points one at a time. This helps to reduce the size of the normal matrix significantly, thereby making the triangulation of a continent-wide scale block feasible and efficient. One of the unique characteristics of our method is the use of a correlation model linking the attitude corrections between images of the same pass, which promotes consistency in the attitude corrections. We evaluated the performance of our triangulation method over the Australian continent using the Australian Geographic Reference Image (AGRI) dataset as a reference. Both a free adjustment, using only the pointing information of the Landsat 8 spacecraft, and a constrained adjustment using the AGRI as external control were performed and the results compared. The Australian block’s horizontal accuracy improved from 15.4 m to 3.6 m with the use of AGRI controls and from 15.4 m to 8.8 m without the use of AGRI controls.


Introduction
With the recent advancements in remote sensing and data processing technologies, there has been a substantial increase in the use of medium-and high-resolution satellite data to monitor and assess the changes in the global landscape. The scientific users of today have access to more data that are at a higher resolution in spatial, spectral, and temporal coverage than ever before. This can be attributed to an increase in the number of remote sensing missions across the globe which make the data available at very low to no cost to the end users. Although the large volume of data helps to analyze and better understand land cover changes over time, it increases the disparity between the data, which necessitates additional data preparation work before valuable information can be extracted. There is a growing trend, however, among imagery providers to generate analysis-ready data (ARD), with a goal of reducing the magnitude of data processing and providing data that are consistent and of the highest scientific standards.
To achieve this goal, the task of georeferencing has now assumed greater importance, and it is essential to ensure that orthorectified products generated over time are well registered to the ground. For Landsat, the geolocation operation is performed by registering to a global reference dataset, referred as the Global Land Survey (GLS), using correlation-based image matching techniques [1,2]. Although the current GLS version, circa 2000 (GLS2000), is a major improvement over earlier versions, large geolocation errors and inconsistencies between adjacent scenes have been observed in some of the GLS data. The GLS2000 in general is accurate to 25 m root mean square (RMS) error, but studies using Landsat 8 with an absolute geolocation accuracy of 18 m at 90% circular error (CE90) suggests that there are several scenes with geolocation errors on the order of 50 to 100 m [3,4]. The large errors observed in the GLS correlates to the regions where accurate ground controls were unavailable during GLS construction. The large inaccuracies and inconsistencies motivated us to work on improving the GLS dataset using better methods, such as space-based triangulation adjustment.
The bundle adjustment techniques using space-based triangulation have been used successfully for validation and have produced consistent and accurate datasets [3,[5][6][7][8][9][10][11][12][13][14]. Some of these methods use traditional bundle adjustment techniques used in aerial photogrammetry by adjusting the interior and exterior orientation parameters ( [11,15]), while others have refined this approach by using rational polynomial coefficients ( [13]) and with additional constraints provided from a digital elevation model (DEM) [8,9]. Many of these studies focused on using bundle adjustment techniques for improving the geoposition of the data but were limited to a few images or over a small geographic area. The study by GeoScience Australia (GA) ( [5]) is an exception, as they produced a georeferenced dataset at a continental-wide scale. The GA produced a georeferenced dataset over Australia by using a long-strip adjustment technique that facilitated a refinement in the georeferencing process of orbit and attitude parameters of Advanced Land Observing Satellite (ALOS) Panchromatic Remote-Sensing Instrument for Stereo Mapping (PRISM) images [5]. This dataset is referred to as the Australian Geographic Reference Image (AGRI). Their studies suggest that the long-strip adjustment technique using a limited number of accurate ground control points can georeference to about 1-pixel accuracy with an expected RMS error close to 2.5 m and a 90% CE of 5.5 m. This is a significant advancement in the process of generating reference imageries using space-based strip adjustment, but implementation of such a technique directly to improve the GLS dataset is not feasible. For example, the strip adjustment technique corrects for bias errors in the attitude and ephemeris, but does not account for the error in the time-varying aspect of these parameters. The bias-only adjustment, though reasonable for the relatively small orbital pass duration over Australia, can introduce large errors in the geoposition for longer passes. This will invalidate the requirement to have reference datasets that are not only accurate but also consistent with the adjacent regions. Other triangulation-based bundle adjustment methods involve a large number of accurate ground control points and use rational polynomial coefficients. This will increase the number of adjusted parameters and the corresponding complexity of the normal matrix especially over large geographical regions. Our approach differs from others by applying bundle adjustment techniques using a physical exterior orientation model to improve the geolocation errors in the triangulated ground points. The purpose of our satellite block triangulation algorithm is to integrate the information from a large number of image point observations, sparse ground control information, and the a priori knowledge of the spacecraft position and attitude to yield a geometrically consistent and minimum variance solution. While most other triangulation applications are directed at correcting the triangulated images, we are using triangulation in a geopositioning application to create an accurate and consistent control point set.
In this paper, we have provided a detailed description of our space-based triangulation algorithm that was developed and implemented to improve the geopositional accuracy of the GLS dataset using Landsat 8 images. The remainder of this paper is organized into three sections. The overview of our algorithm, theoretical construct and structure of the normal matrix used in the triangulation procedure, and the method used to reduce and solve the large normal matrix using sequential elimination techniques are introduced and discussed in Section 2. The application of this algorithm to improve the GLS over Australia is discussed in Section 3, and the overall conclusions from our study are given in Section 4.

Methods
The bundle adjustment using triangulation-based methods involves iteratively adjusting the spacecraft position, velocity, attitude, attitude rate corrections, and ground point coordinates so that the residual errors in the measured image points are minimized. The observation equation or look-point equation, which relates the ephemeris, attitude, and ground control point data to the image observations, is non-linear. Our approach uses a Taylor series approximation with numerically computed partial derivatives to linearize the look-point equation. The linearized equation is solved iteratively using an initial approximation, a priori covariance information, and computed partial derivatives. Each iteration yields a set of incremental correction parameters that are accumulated over all iterations. As the iterative solution converges, the incremental adjustments computed in each iteration grow progressively smaller. The iterations are halted when the adjustments to the parameters are smaller than a specified convergence threshold. In essence, the method solves a non-linear least squares problem using a priori knowledge of the spacecraft and ground control points to minimize the residual errors from image-based observations of the ground control points. The mathematical formulation and description of our method is described in the following sections.
A pictorial representation of the space-based triangulation procedure is shown in Figure 1. The orbital tracks prior to the triangulation for the two passes, pass 1 and pass 2, are shown in solid red and blue curves, respectively. The corresponding dotted lines (red, blue) denote the orbital track after the triangulation procedure. The solid red and blue line of sight vectors from the spacecraft to the ground represent the viewing geometry of a specific ground point observed in the image, prior to the triangulation process. The adjusted line of sight vectors to the ground point after the triangulation are denoted by dotted lines. The pre-triangulated positions of the ground points are shown in black dots, and their corresponding locations estimated from the spacecraft's pre-adjustment attitude and ephemeris are shown by red and blue dots for pass 1 and pass 2, respectively. The same points in image 2 are shown as red and blue circles for passes 1 and 2, respectively. The green stars represent the triangulated positions of the ground points, and the black arrows indicate the positional offsets (∆G) for the ground points after triangulation. In general, a spacecraft flying in an orbit (solid red curve) images a ground point (black dot) at a location offset from the ground point location (geodetic offsets), due to the inaccuracies (random errors) in the spacecraft's position and orientation. In cases where the ground point locations are also inaccurate, the true position of the ground points cannot be determined accurately due to the uncertainties in the observations. The space-based triangulation process reduces the uncertainty in the estimation of the true position of the ground point by using additional information provided by the other passes or images. In this case, pass 2 observes the same ground point but locates it on the ground (blue dot) at a different place than pass 1. Using the a priori information of the spacecraft's ephemeris, attitude, and the ground control point, the minimum-variance-based bundle adjustment technique simultaneously solves for the new position of the ground control point that minimizes the overall error in the spacecraft and ground parameters. In Figure 1, the corrections to the spacecraft parameters are denoted by (∆P 1 , ∆V 1 , ∆α 11 , ∆ᾱ 11 ) and (∆P 2 , ∆V 2 , ∆α 21 , ∆ᾱ 21 ) for passes 1 and 2, respectively, and the corresponding adjustments to the ground parameters are denoted by ∆G. The advantage of the bundle adjustment is that the ground points that are observed only in a single image (and pass) are also improved by adjusting the spacecraft's position and orientation parameters. Figure 1. A pictorial representation of the space-based triangulation procedure. The orbital tracks prior to the triangulation for the two passes, pass 1 and pass 2, are shown in solid red and blue curves, respectively. The corresponding dotted lines (red, blue) denote the orbital track after the triangulation procedure. The solid red and blue line-of-sight vectors from the spacecraft to the ground represent the viewing geometry of a specific ground point observed in the image, prior to the triangulation process. The adjusted line-of-sight vectors to the ground point after the triangulation are denoted by dotted lines. The pre-triangulated positions of the ground points are shown in black dots, and their corresponding locations estimated from the spacecraft's pre-adjustment attitude and ephemeris are shown by red and blue dots for pass 1 and pass 2, respectively. The same points in image 2 are shown as red and blue circles for passes 1 and 2, respectively. The green stars represent the triangulated positions of the ground points, and the black arrows indicate the positional offsets (∆G) for the ground points after triangulation.

Data Synthesis
This section discusses the formulation and the procedure used to construct and extract the necessary information/data to set up the observation equations for least squares.

Solution Parameters
The solution parameters or unknowns in our bundle adjustment process consist of 6 ephemeris correction parameters (3 position and 3 velocity) for each satellite pass, 6 attitude correction parameters (3 bias and 3 rate) for each image, and 3 position correction parameters for each ground point measured as shown in Equation (1). In this context, an image refers to a single Worldwide Reference System (WRS-2) Landsat scene, while a pass consists of an entire descending orbit and may contain one or more images, which need not be contiguous. The WRS-2 is a global notation system used to reference and catalog Landsat data. It enables the user to identify a specific geographic region of the world in the Landsat archive by the corresponding WRS-2 designation, which uses two numbers namely, path and row. The path numbers identify the nominal satellite orbital tracks, while the row numbers relate to the latitudinal center line of a frame of imagery or scene. The "pass" in this paper refers to an image acquisition interval by Landsat over a WRS path on a particular day, while an image refers to a WRS scene within a pass, which is approximately 24 seconds of contiguous data. The pass corrections are limited to one set of ephemeris corrections to account for the fact that the ephemeris errors are highly time correlated due to the dynamics of the spacecraft orbit. Furthermore, restricting the ephemeris correction parameters to a single set per pass also helps to decorrelate the ephemeris corrections from the attitude corrections. Attitude corrections are introduced for each image to account for the more rapidly varying attitude knowledge errors. But the parameters for all images from the same pass are linked using correlation observations to provide along-track continuity in the attitude correction solution. This unique approach of linking attitude corrections within a pass using a time-based correlation model limits both the magnitude and direction of attitude correction of each image. The attitude correlation model is explained in Section 2.2.2. Each image has an associated scene center time (tsc j ), and each pass has an associated pass center time (tpass j ), which is calculated as the median of the scene center times (tsc j ) for the scenes in the pass. Note that there will be one or more image observations associated with each ground point, which may come from adjacent images of a single pass and/or images from adjacent passes.
where ∆P k and ∆V k are the ephemeris position and velocity corrections, respectively (one set per pass "k") in the along-track (AT), across-track (XT), and radial (Radial) directions; ∆α k and ∆ᾱ k are the attitude bias and rate corrections, respectively (one set per image "j") in the roll, pitch, and yaw orientations; ∆X g n refers to the ground control point corrections (one set per ground point "n") in the Earth-centered-Earth-fixed (ECEF) X, Y, and Z coordinate axes. The correction terms as defined in Equation (1) are initialized to their a priori value of zero. These a priori parameters have associated a priori covariance information, which is provided by the user as parameters of the triangulation solution in the form of a priori weights on the ephemeris, attitude, and ground control data. The covariance information for the ephemeris and attitude is dependent on the characteristics and accuracy of the spacecraft. For example, the Landsat 8 spacecraft uses an on-board global positioning system (GPS) to estimate its position and velocity, and has sophisticated and advanced gyro systems with low drift errors, which are more accurate and precise than the systems used in earlier Landsat missions (Landsat 5, Landsat 7). Therefore, it is important to use the appropriate a priori covariance information that is tailored to the specific satellite system. Similar to spacecraft position and orientation, the ground point covariance data can also be varied depending on the accuracy of the data source and the type of ground point. For example, GPS surveyed points can be assigned a smaller a priori variance (larger a priori weight) than control points transferred from orthorectified images. Similarly, within the triangulation procedure, control points may be assigned a smaller a priori variance than a tie point. As in the case with spacecraft covariance information, the covariance for ground points should be chosen based on the expected accuracy of the control point data. An underestimate or overestimate of the a priori covariance can lead to an increased number of iterations within the least squares and also produce larger residuals, resulting in poor correction estimates for the pass, image, and ground points. In this study, we used the existing GLS ground control points used by the Landsat product generation system [3,16,17] as control points with sufficiently large a priori variance as they are expected to be less accurate. The a priori weights used in this study are discussed in Section 3.

Image Point Observations
In this study, we used the Image Assessment System (IAS), a tool developed by U.S. Geological Survey (USGS) for Landsat calibration [17,18], to measure the location of the ground control points in the images. The control points are measured in Landsat 8 images that have been geometrically corrected with compensation for terrain relief effects, but without reference to ground control-a systematic terrain-corrected (L1GT) image product. The resampling grid used to generate the L1GT image allows us to map the measured control point locations back to the original unresampled sensor space. This allows us to recover the time of observation and the corresponding spacecraft position and velocity. The measured coordinates in the Earth-referenced L1GT image and the corresponding height from the DEM used to create the L1GT image provide the information needed to construct the apparent ECEF ground point position. We construct the image observation vector using the spacecraft position at the time of observation and the apparent ground point position, rather than explicitly reconstructing the sensor viewing geometry from the sensor model. For each observation point i, several data elements, as shown in Table 1, are collected to formulate the coordinate axes, the observation vector, and to compute the observation covariance as shown in Equation (2)-(4). For each observation, we record the observation vector (x i ), time of observation (tobs i ), spacecraft position and velocity (P i , V i , respectively), and observation covariance (Σ i ). The orbital coordinate system axes are given as Using this, we can construct the observation vector as Each image measurement adds two observations, one in the along-track and the other in the across-track directions. Ground points with multiple image observations from different paths can be positioned in the 3 dimensions (horizontal and vertical), whereas the points with single image observations rely on the a priori knowledge of the height to update the horizontal position. The observation covariance is estimated numerically by perturbing the line i and samp i (by adding 1) without changing the tobs i . For example, x i is calculated for (line i + 1, samp i ), and the partial of x i with respect to line is calculated by taking the difference of x i from x i . Similarly, the partial derivative of x i with respect to sample is calculated by perturbing samp i .

Ground Point Information
All the ground control points used in the study are specified in WGS84 ellipsoidal latitude, longitude, and height. The ground point position vector in the Earth-centered-Earth fixed coordinate system (ECEF), X g n , is computed from the geodetic latitude (φ), geodetic longitude (λ), and height (h) using the standard transformation [18], as shown in Equation (5).
The ground point position vector is given as where N = a √ 1−e 2 sin 2 φ , and a and e are the WGS84 ellipsoidal semi-major axis and eccentricity, respectively. The covariance for the ground point vector in ECEF (Σg n ) can be computed from the geodetic covariance matrix (Σφ) as where σ φ is the latitude standard deviation based on the accuracy of the ground control point (GCP); σ λ is the longitude standard deviation based on the accuracy of the GCP; σ h is the height standard deviation based on the accuracy of the GCP.
The Jacobian matrix (W) is determined directly by taking the partial derivative of the trigonometric functions in Equation (5), as shown in Equation (7).
For each ground control point, we record the ground control position vector (X g n ) and its covariance matrix (Σg n ).

Algorithm Overview
As discussed earlier, the space-based bundle block adjustment technique linearizes the non-linear problem using a Taylor series expansion to solve for the correction estimates in attitude, ephemeris, and ground control points using a least squares approach. The least squares problem is formulated as Nδ = C, where N is the normal equation matrix, δ is the vector of unknowns as discussed in Section 2.1.1, and C is the right-hand-side vector. Note that in an iterative solution such as this, the unknowns are the current iteration's incremental adjustments to the net correction parameters. These incremental corrections should be driven to zero (or convergence threshold) as the iterative solution converges. At the beginning of each triangulation solution iteration, the normal equations are initialized with a priori "observations" that capture our knowledge of the accuracy of the ephemeris and attitude data. These observations constrain the magnitude of the ephemeris and attitude corrections and also provide the linkages between attitude correction parameters from the same pass. Once the normal equations are initialized, the ground points are processed one at a time. All image observations for a given ground point are processed, and that ground point is eliminated (temporarily) from the triangulation solution before proceeding to the next ground point. Once all the ground points have been processed, the "reduced" normal equations, containing only the ephemeris and attitude correction unknowns, are solved, and the resulting corrections are used in a back-substitution procedure to calculate the (previously eliminated) ground point correction unknowns. The advantage of this parameter reduction/elimination technique is that the size of the normal equation matrix will be considerably smaller, and therefore, its inverse can be computed with limited memory requirements in a reasonable time. The formulation of this ground point elimination procedure is discussed in detail in Section 2.2.3.
We begin the least squares process by gathering all the observations for each ground control point, which are then used to construct the normal equation matrix as described below.
For each ground point, n: Update the initial ground position with the current correction estimate: X g n = X g n + ∆X g n .
(2) For each image point i, viewing ground point n, in image j and pass k: (a) Compute the local ECEF to orbital coordinate system rotation matrix, M i , using the spacecraft state vector for this tie point, P i , V i (see Equation (2)).
(b) Update the spacecraft position and attitude with the current corrections. Note that the corrections ∆P k and ∆V k are in the orbital coordinate system; therefore, they must be converted to ECEF to be applied to the initial state vector.
Compute the updated ECEF to orbital ( M) and attitude correction ( T) matrices: The attitude correction rotation matix (T) for roll (r), pitch (p), and yaw (y) is given as (d) Construct the predicted three-dimensional observation vector d(u i , ν i ), which expresses the image measurements, u i , and observation residuals, ν i , as a function of the point, pass, and image unknowns for observation i.
We can eliminate the unknown scale factor "s" by dividing the first and second element by the third element to construct a two-dimensional observation vector x i Construct the line of sight/look-point equation by combining the measured and predicted observation vectors.
where ν i is the observation residual vector for observation i.
Linearizing the function (Equation (11)) via a Taylor series expansion gives the following observation equation: where ν i represents the linearized observation residual vector, and A represents the Jacobian matrix or matrices of partial derivatives with respect to the ephemeris corrections (A 1ik ), attitude corrections (A 2ij ), and ground point corrections (A 3in ) for image point i, ground point n, in image j and pass k, as shown in Equation (13).
These partial derivatives are computed numerically by perturbing the correction parameters one at a time and computing the resulting effect on the predicted observation vector, for example: (g) Construct the normal matrix: Assuming the normal equation Nδ = C, the normal matrix N can be partitioned as shown in Equation (15), where δ 1 , δ 2 , and δ 3 correspond to the ephemeris corrections, attitude corrections, and ground point corrections, respectively. The structure of this matrix formulation is discussed in Section 2.2.1.
For the ground point n, image point i, in image j, and pass k, the contribution of this observation to the normal equations can be estimated from the Jacobian matrix (see Equation (13)) as shown below in Equation (16), where Σ i is the image-space covariance of the observation (see Equation (4)). The numbers in the parentheses indicate the dimensions of the matrix. Note that this accounts only for the contribution of the specific observation and does not include the a priori contributions.
The individual observation contributions to the normal matrix are combined with the a priori contributions, which are discussed in Section 2.2.1. In the traditional least squares approach, the entire normal matrix is inverted to estimate the parameter adjustments/correction in each iteration until convergence is achieved. However, as discussed in the next section (Section 2.2.1), the size of the N matrix can quickly become large, even for a small geographic region, and cannot be easily inverted. We mitigated this issue by eliminating the ground points one at a time, the details of which are discussed in Section 2.2.3.

Structure of the Normal Equations
If we have K passes, J images, and N ground points in the bundle adjustment solution, then there will be a total of 6K + 6J + 3N unknowns, giving the matrices in the normal equation Nδ = C the following dimensions: As there should be multiple ground points per image, a typical solution will have: N > J > K, which tends to make the number of ground points drive the size of the normal equations to be solved. Fortunately, as explained below, it is not necessary to solve the full (6K + 6J + 3N) by (6K + 6J + 3N) problem simultaneously. The ground point unknowns can be added to and eliminated from the solution sequentially, making it possible to reduce the dimension of the simultaneous solution to (6K + 6J) by (6K + 6J). The ground point unknowns can be solved subsequently by back-substitution. To see how this is done, consider the structure of the normal equations as shown in Equation (15). The elements labeled "1" correspond to the ephemeris corrections (one set of 6 for each pass k), the elements labeled "2" correspond to the attitude corrections (one set of 6 for each image j), and the elements labeled "3" correspond to the ground point corrections (one set of 3 for each ground point n). The diagonal elements of N, N 1 and N 3 , are block diagonal, while N 2 is block semi-diagonal, as shown in Equation (18). N 1 has one 6-by-6 non-zero diagonal component for each pass, and conceptually, N 3 has one 3-by-3 diagonal component for each ground point, although we never form the complete N 3 matrix. N 2 has one 6-by-6 diagonal component for each image and a 6-by-6 off-diagonal block linking each pair of images from the same pass.
The matrix structure above assumes that images j and m are from the same pass. The N 1ik , N 2ij , and N 3in terms are the contributions from image observation i as shown in Equation (16). The N 10 , N 20 , and N 30 terms are the contributions from the a priori weights applied to the ephemeris, attitude, and ground point corrections, respectively. As such, they are equal to the inverses of the a priori correction term covariance matrices, as shown in Equation (19) below. Note that N 10 and N 20 are the same for all passes and images, respectively, whereas N 30n is indexed by ground point number (n), allowing the ground points to be individually weighted. The N 2j and N 2jm terms are used to link the attitude corrections for images from the same pass by applying attitude correlation between two image pairs, the details of which are discussed in Section 2.2.2.
The off-diagonal matrices of the normal matrix N in Equation (15) are all sparse matrices. The M matrix with dimension 6K by 6J has 6 by 6 non-zero components along the row corresponding to pass k, only in those columns that correspond to images j, that are part of pass k. The Q matrix with dimension 6K by 3N has 6 by 3 non-zero components along the row corresponding to pass k, only in those columns that correspond to ground points n, that are imaged in pass k. Similarly, the R matrix with dimension 6J by 3N has 6 by 3 non-zero components along the row corresponding to image j, only in those columns that correspond to ground points n, that are viewed in image j.
The M ikj , Q ikn , and R ijn terms are the contributions from image observation i, as shown in Equation (16). The δ 1 vector (6K by 1) contains the ephemeris corrections, the δ 2 vector (6J by 1) contains the attitude corrections, and the δ 3 vector (3N by 1) contains the ground point corrections for a given iteration. The C vectors contain the right-hand-side values for the ephemeris, attitude, and ground point corrections. The C 1 , C 2 , and C 3 matrices have the same dimensions as δ 1 , δ 2 , and δ 3 .
The C 1ik , C 2ij , and C 3in terms in the above equations are the contributions from image observation i as defined in Equation (16). The C 10k , C 20j , and C 30n terms in Equation (22) correspond to the a priori weights on the corrections (i.e., observations that the corrections are equal to zero). Note that unlike the corresponding a priori contributions to N, the C terms are functions of the current (net) correction values so they are indexed by pass, image, and ground point. The C 2j is the contribution for the scene-to-scene correlation observations defined in Section 2.2.2.

Attitude Correction Correlation Observations
Landsat images separated by a short duration within a pass are highly correlated in their attitude knowledge errors. Therefore, in the triangulation solution, although each image has separate attitude correction parameters, the correlation between the attitude observations should be established. This is achieved by building an attitude correlation model to estimate the image covariance for image pairs from the same pass. Each image pair in the same pass yields a symmetric pair of observations: (a) the correction for image j predicted from the correction for image m and (b) the correction for image m predicted from the correction for image j. These correlation observations take the following form, as shown in Equation (23), where I refers to the identity matrix.
Note that the signs on the terms containing ∆t in B 2 and D jm are inverted as compared to the corresponding terms in B 3 and D mj because we have maintained a consistent definition of ∆t in all of these equations. The correlation observations are weighted based on the time difference between the scene centers. The weight matrices (W jm and W mj ) are estimated by assuming the attitude correction parameters to be realizations of a random process, a(t), the attitude state vector, with zero mean and a user-provided a priori covariance following a Gaussian distribution, i.e., a(t)∼Gaussian(0, Σ). We further assume that a(t) is a Gauss-Markov process, which implies that the covariance of the two realizations separated in time is given as shown in Equation (24).

E[a(t)a T (t + ∆t)]
and τ is a user-provided correlation time constant parameter.
Based on the definition of the attitude correction model, a(t) evolves with time as a(t + ∆t) = Φ(∆t)a(t) + w(t), where Φ(∆t) = I ∆tI 0 I , and w(t) is process noise.
Our correlation observations link the states associated with each pair of images in a pass by predicting a(t + ∆t) from a(t) using this propagation model. The covariance statistic of this prediction is what we are interested in, and that is given in Equation (26).
The weight matrix W jm , shown in Equation (27), is the inverse of the above matrix (Σ ∆t ).
Using the weight matrices, contributions of the correlation observations to the normal equations are computed once for each iteration by analyzing the set of images, j = 1, ..., J k , in each pass k. For each image pair (j, m), the contributions to the normal equations are shown in Equation (28). Note that each image pair (j, m) yields contributions to the diagonal blocks corresponding to images j and m, and that the individual contributions from each image pair are summed, whereas each image pair (j, m) also yields a unique off-diagonal block at row j, column m, in the block-structured N 2 matrix (as well as its transpose at row m, column j). The correlation observations also generate contributions to the C 2 vector, which are shown in Equation (29). As in the case for the contributions to the diagonal blocks in the N 2 matrix, the individual contributions to C 2 from each image pair (j, m) are summed.

Eliminating Ground Point Unknowns
As stated earlier, we process and eliminate one ground control point at a time to reduce the size of the normal matrix. Once all of the image observations i corresponding to ground point n have been added to the normal equations, point n can be eliminated from the solution as discussed below. (1) First compute the inverse of the 3 × 3 N 3n matrix: N −1 3n . (2) Compute and save (for future back-substitution): For each pass that viewed point n, k = 1, ..., K n : Subtract this term from the row k block of C 1 . (c) Compute and save (for future back-substitution): For each pass that viewed point n, m = 1, ..., K n : Subtract this term from the row m, column k block of N 1 .
For example, if only two passes, k and m, view a point, the N 1 and C 1 matrices would be adjusted as Similarly, for each image that viewed point n, j = 1, ..., J n : Subtract this term from the row j block of C 2 . (c) Compute and save (for future back substitution): For each image that viewed point n, m = 1, ..., J n : Subtract this term from the row m, column j block of N 2 .
(e) For each pass that viewed point n, k = 1, ..., K n : Subtract this term from the off-diagonal block matrix M kj . (iii) Note that for passes k other than the one containing image j, the contents of M kj will initially be zero, but this will change as the process of ground point elimination progresses.
The computations as shown above eliminate ground point n from the solution. The ground points are sequentially processed and immediately eliminated to create the reduced normal equations, as shown below in Equation (31). The ground point corrections can be determined using the back-substitution method as shown in Equation (32).
Having solved for all the unknowns for this iteration, the net corrections are updated as in Equation (33). The iterations are continued until the incremental corrections are sufficiently small.

Results and Discussion
The USGS Landsat processing system uses the GLS2000 dataset as a geometric reference and extracts ground control points (GCPs) in the form of image chips (64 × 64 pixels) centered at well-defined points of interest [16]. These image chips, or GCPs, are used to register and generate Landsat terrain-corrected products. In general, the GLS2000 dataset is accurate globally to within a 25 m RMS error, but there are regions in certain parts of the globe where the absolute positional error can exceed 40 m. This is especially true in places where accurate ground control positions were unavailable when the GLS was created [3]. With the help of the satellite block triangulation technique discussed in this work, we decided to improve the absolute accuracy of the GCPs derived from the GLS2000 dataset using the improved pointing accuracy and performance of the Landsat 8 spacecraft.
In this work, we applied our triangulation method to perform a continental-level block triangulation over Australia. The Australia block contains all path/row locations in Australia spanning WRS paths 088 through 115 and rows 068 through 090, a total of 394 path/row scenes. Although none of the GLS2000 scenes in this area exhibited geodetic offsets (absolute positional accuracy) of over 50 m, there were several places where offsets greater than 30 m were observed. Geoscience Australia (GA) performed a comparison of the GLS dataset with their highly accurate Australian Geographic Reference Image (AGRI) [5]. Their assessment also identified a few regions where the discrepancies between the GLS and AGRI exceeded 30 m. To improve the consistency between the GLS and AGRI framework, we performed a triangulation-based adjustment of the GLS GCPs.
We evaluated the performance of our triangulation-based bundle adjustment technique for two cases: (a) unconstrained or free adjustment, and (b) AGRI-constrained adjustment. In the free adjustment case, the absolute pointing information from the Landsat 8 Operational Land Imager (OLI) scenes was used to improve the position of the GLS GCPs. This provides an unconstrained triangulation solution in the sense that the Landsat 8 pointing accuracy is used to improve the GLS GCPs without using any external source of ground control. The GLS control points were assumed to have poor positional accuracy and were given very low weights. Table 2 shows the a priori standard deviations (inversely related to the weights) used in the triangulation process. In the case of the AGRI-constrained adjustment, the GCPs extracted from the AGRI reference images were used as control points (held fixed by using high weights), and the constrained bundle adjustment provided the updated positions for the GLS GCPs, which, as in the free adjustment case, were given much lower a priori weights. For both cases, the ground point measurements in the OLI images were performed using the Image Assessment System (IAS) tool [18]. As shown in Table 2, the horizontal standard deviations for the GLS ground points were set to 10,000 m for both cases. The GLS GCP network was quite dense, with hundreds of points per Landsat WRS-2 scene. This dense GCP coverage is used to ensure that even predominantly cloudy images can be registered. A typical scene will generate about 1000+ observations (correlates 1000+ GCPs). As a result, the effective variation allowed by the GCPs as a group is reduced by √ n, i.e., the net a priori GLS GCP standard deviation for a typical scene will be 10,000 √ ≈ 300 m. Although 10,000 m does not reflect the accuracy of an individual GLS GCP, it was necessary to use a very low a priori GCP weight to prevent the GCPs in aggregate from overly influencing the solution, given how numerous they were. In the case of the AGRI-constrained adjustments, the AGRI control points were allowed to vary by no more than 5 m, which was the expected accuracy of the AGRI dataset [5]. Since the AGRI provides only horizontal positions, we used Shuttle Radar Topographic Mission (SRTM) data as the elevation reference. These data were in general accurate to better than 30 m, and so the standard deviation for the GLS ground point height was set to 30 m [19,20]. The highly accurate GPS sensor onboard the Landsat 8 (L8) spacecraft enabled us to provide a small standard deviation of 5 m for the ephemeris and 1 mm per sec for the velocity. Similarly, the accurate attitude sensors on the L8 allowed us to use 10 micro-radians (1σ) as a measure of allowable corrections to the attitude for each image. The standard deviation estimates for the AGRI GCPs, ephemeris, and attitude were chosen to be small so that the corrections were predominantly in the GLS GCPs, thereby transferring the accuracy of the OLI scenes (free adjustment) and the AGRI control points (AGRI constrained) to the GLS. Furthermore, our 1σ estimates for the ephemeris and attitude parameters were consistent with the established pointing accuracy of L8 based upon on-orbit characterization, which indicated an absolute accuracy of 18 m CE90 [4]. The GCPs were identified in the image using normalized grayscale correlation techniques. The expected uncertainty in this image measurement was used as a measure of standard deviation for the observation point.Typically, well-defined GCPs are located to better than quarter of a pixel, and therefore, we chose 5 m as the 1σ estimate to account for correlation-related errors in the triangulation solution [4]. Table 2. Standard deviations used to determine the weights in the free adjustment and Australian Geographic Reference Image (AGRI)-constrained triangulation solutions. The AGRI data provided by GA consists of 2.5 m ground sample distance (GSD) panchromatic images from the ALOS PRISM sensor. These images were downsampled to 15 m GSD and mosaicked to generate reference images that cover the geographic extent of Landsat 8 (WRS-2) scenes. The anchor sites were used as control sites in the development of the GLS2000 dataset and provide a good distribution across Australia. The validation sites were chosen to include any regions that were identified by GA to have large errors, and additional sites were added for a more complete geographic distribution. In this study, GCPs (image chips) extracted from the AGRI reference images over the anchor sites were used as AGRI ground control points in the triangulation process. Triangulation of the Australia block included 394 images (scenes) from 266 passes (intervals) to yield updated positions for 134,745 GLS GCPs. For the AGRI-constrained case, an additional 10,897 AGRI control points were also part of the triangulation, but these points were held as a fixed control to constrain the triangulation solution. Due to the large size of the Australia block, only a single scene from each path/row was used in the triangulation. As mentioned in the methodology section, the Landsat images within a pass were highly correlated in their attitude estimates over short duration and were linked using the attitude correlation model described in Section 2.2.2. From our independent study using different trials, we observed that the attitude errors typically decorrelate after 60 s, and therefore, the attitude correlation time constant parameter (τ) was set to 60 s. Table 3 shows the number of observations, ground control points, images, and passes in the triangulation solution for the two cases. In both cases, the triangulation method converged to a solution in 2 iterations and reduced the overall observation error by a factor of 10. Improvements to the GLS GCPs' positions were assessed using two validation methods provided by the IAS: (a) the geodetic accuracy assessment procedure, and (b) the image-to-image (I2I) characterization procedure. In the first method, we estimated the geodetic offsets between the updated GLS GCPs and the Landsat 8 scenes that were processed using only spacecraft pointing and relief adjustments (systematic terrain-corrected product). Details of this IAS characterization procedure can be found in the Landsat 8 documentation [18]. For this validation, we used the same Landsat 8 scenes that were used in the triangulation process to generate the GLS GCP observations. In the second method, orthorectified Landsat 8 products generated using the updated GLS GCPs (as control points) were compared directly against the AGRI reference images at the control and validation sites using the I2I characterization procedure in IAS [18]. Unlike the first method, this comparison directly provided an estimate on the absolute accuracy of the GLS GCPs in comparison with the AGRI dataset.

Items/Description
Geodetic offsets were measured for the original unadjusted GLS control, the control resulting from the free adjustment, and the control from the AGRI-constrained adjustment. Summary statistics for the geodetic offsets from these three cases are shown in Table 4. The geodetic offsets of the Landsat 8 systematic terrain-corrected products with respect to the GLS control prior to the triangulation procedure are shown in Figure 2, where the anchor sites are shown in black outline. A smaller geodetic offset indicates that the GLS GCPs were consistent with the Landsat 8 scenes' positional accuracy, whereas a larger number suggests inconsistencies between the two. The main advantage of this plot is that it helped identify regions where the inconsistencies were larger than 1 multispectral pixel (30 m). GCPs in these regions (shown as white ellipses) were highly likely to have poor positional accuracy and should be improved after the triangulation procedure. The vector plot on the right shows the magnitude and direction of the geodetic offsets (residuals) for each of the GCPs for a particular path/row (90/82) whose net geodetic offset was greater than 30 m (red). The magnitude of the residuals was scaled by a factor of 300 for visual purposes. In this plot, the downward-pointing arrow indicate a large systematic bias in the along-track direction. The mean geodetic offsets for the scene (Landsat scene id: LC80900822013188LGN01) used in the triangulation procedure were found to be -3.1 m and -31.9 m in the across-and along-track directions, respectively. The geodetic offsets after the triangulation procedure are shown in Figures 3 and 4 for the free and AGRI-constrained adjustment cases, respectively. Comparing the pre-and post-triangulation geodetic results, it was observed that the regions identified with poor geometric accuracy (white ellipses) in the original GLS had improved significantly in both the free and AGRI-constrained cases. Furthermore, offsets for many of the scenes in the block were within 6 m for the free adjustment case. This was to be expected as the free triangulation adjustment process essentially improves the entire block using the pointing knowledge of these scenes. In the AGRI-constrained case, the offsets for most of the scenes were between 6 m and 12 m, which was consistent with the expected accuracy of the Landsat 8 pointing knowledge and AGRI reference. For one of the scenes in the block (near the center of the block), the offset observed was larger than expected (24-30 m) for the free adjustment case, but the same scene was also found to show large offsets in the AGRI-constrained case (see Figure 4). This suggests that the error in this path/row was mainly due to poor pointing knowledge of the specific L8 scene, and not in the improved GCPs. This clearly demonstrates that the triangulation procedure has the ability to make an entire block consistent without introducing scene-dependent bias into the final adjustments.  The second validation method used the I2I characterization procedure to assess the improvements from the triangulation technique. Figures 5-7 show the geometric error between the AGRI reference and the Landsat 8 orthorectified scenes (using GLS GCPs) prior to the triangulation, after free-adjustment triangulation, and after AGRI-constrained adjustment, respectively. The summary statistics of this validation method for the original GLS control and for the two triangulation solutions are shown in Table 4. In Figure 5, it was observed that the regions with large geometric offsets (>30 m) are the regions where the geodetic characterization procedure also showed large offsets (see Figure 2). Similar to the geodetic offsets, large offsets in I2I indicated regions where the GLS GCPs should be improved. Unlike geodetic offsets, which are dependent on the pointing accuracy of L8 (18 m), the I2I offsets can directly provide the absolute error of the GLS GCPs by comparing with the more accurate AGRI reference (5 m). Comparing the pre-triangulation adjustment with the free adjustment results (see Figures 5 and 6, Table 4), it was evident that the regions with large residuals were improved significantly and that the majority of the sites were within 12 m, even without using any ground controls in the block. The net RMS for the 37 anchor sites improved from 15.4 m in the original GLS to under 9 m after free adjustment, and the largest residuals in the block reduced from 36.8 m to 13.1 m. The RMS and the maximum offset in the block for the anchor and validation sites were similar, which indicated that the original GLS triangulation effectively distributed the control accuracy from the anchor sites throughout the Australia block. In the validation for AGRI-constrained triangulation, the RMS for the anchor sites improved from 15.4 m to 3.1 m, and the largest offset in the block reduced from 36.8 m to 5.7 m. As in the free adjustment validation, the RMS was consistent for both the anchor and validation sites, and the overall offsets for all the validation scenes except for one were under 6 m. A limited number of external ground control points (37 AGRI scenes versus 394 scenes in the block) were enough to adjust the entire block without introducing any geometric artifacts in the block, as evident in Figure 7. The validation results presented in this section have clearly demonstrated that our triangulation method can adjust a continent-level block seamlessly and consistently both with and without the need for additional ground control to improve the block. Although there was a significant improvement with the use of ground control points, our triangulation method can make use of the accurate pointing knowledge of the spacecraft alone to generate a consistent, global geometric reference GCP set.  Residual offsets between Landsat 8 orthorectified products using GLS GCPs after free triangulation adjustment and the AGRI reference data using the image-to-image (I2I) characterization procedure. The anchor sites are indicated by black outline. Note that by using the free adjustment triangulation, most of the sites are within 12 m with respect to the AGRI reference even without using the AGRI control in the triangulation. Residual offsets between Landsat 8 orthorectified products using GLS GCPs after AGRI-constrained triangulation adjustment and the AGRI reference data using the image-to-image (I2I) characterization procedure. The anchor sites are indicated by black outline. Note that using the AGRI-constrained solution, almost all of the test sites are within 6 m with respect to the AGRI controls.

Conclusions
In this paper, we have described the space-triangulation-based bundle adjustment algorithm and demonstrated how the algorithm can be used to improve the local and regional inconsistencies in the GLS2000 reference dataset and create a refined, consistent, and accurate control base. In general, the geometric improvements in the positional accuracy and consistency of the data are dependent on the type of sources used in the triangulation. In this study, we leveraged the pointing knowledge of the Landsat 8 spacecraft to transfer its accuracy to the GLS2000 ground control points. One of the unique characteristics of this algorithm is the use of a correlation model linking the attitude corrections between images of the same pass. This promoted consistency in the attitude corrections, which was important especially when data were acquired over long orbital passes. For example, the L8 satellite passing over Siberia will image the ocean for several minutes before it passes over Australia. The attitude will be highly correlated over Siberia as the images are acquired within a short duration from each other. However, the attitude estimates over Australia will be completely decorrelated with the attitude estimates over Siberia. Therefore, it is important to model this correlation over shorter durations to ensure consistency and to model the decorrelation over long durations to suppress the propagation of attitude errors in the bundle adjustment. Our algorithm assigned a priori covariance for all the parameters in the bundle adjustment, which helped to adjust the block without the requirement of external ground controls. This was important, especially when improvements in the dataset were needed over regions where accurate ground control points were unavailable.
The theoretical basis and methods explained in great detail in this paper were implemented to demonstrate the capability to improving the GLS2000 dataset. In this regard, we processed and analyzed the entire Australian continent both with and without the use of external ground controls to evaluate the improvement and consistency of the entire block. We validated the performance of our method using the AGRI reference data, which were accurate to within 5 m CE90. Comparison of the improved GLS2000 using our bundle adjustment method with the AGRI reference for the free solution (without control) indicated a significant improvement in the absolute positional accuracy and internal consistency of the dataset. The relative horizontal accuracy improved from 15.4 m for the block to 8.8 m without the use of external controls. Furthermore, most of the validation sites were found to be within 12 m with respect to the AGRI dataset, which was in good agreement with the expected pointing accuracy of the Landsat 8 data. When the block was adjusted by including the AGRI control points, the horizontal accuracy improved from 15.4 m to 3.6 m. The residuals for all the validation scenes, except one, were within 6 m, indicating block-level refinements both in accuracy and consistency.
In the future, we plan to extend this work to cover all the continental blocks in an effort to generate a global reference dataset that will be used to produce Landsat orthorectified products. To improve the consistency of the Landsat terrain-corrected products with products from other sensors such as Sentinel 2, we plan to tie our bundle adjustment to the Sentinel 2 Global Reference Image (GRI) dataset. The GRI dataset is derived from orthorectified Sentinel 2 cloud-free images and is planned to be used as reference in the Sentinel 2 product generation system [21]. We believe that the Landsat 8 terrain-corrected data, when processed using the improved reference dataset, will align to sub-pixel precision with the Sentinel 2 terrain-corrected data that are processed using GRI. Based on the results from the Australian block using the AGRI controls, we anticipate the Landsat 8 to Sentinel 2 registration accuracy to be on the order of 10 m (2σ) or better once the archived Landsat 8 and Sentinel 2 data have been reprocessed using their respective reference datasets. The USGS has plans to reprocess the entire Landsat archive using the improved Landsat global reference dataset as part of Landsat collection processing (Collection 2). This reprocessing is expected to be complete in early 2020. The schedule for reprocessing the Sentinel 2 archive using the GRI has not yet been released officially. The harmonized global references will help in generating terrain-corrected products that are registered to a common geometric reference. This will allow the remote sensing community to use the products from either sensor without a need for registration or double resampling for their time-series analysis [22][23][24], thereby bringing the scientific community one step closer to the realization of analysis-ready data.