ORTHORECTIFICATION OF WORLD VIEW 2 STEREO PAIR USING A NEW RIGOROUS ORIENTATION MODEL

The use of very high resolution satellite images in photogrammetric works is increasing continuously, especially in works such as orthoimages production or Digital Elevation Models extraction. A large number of new photogrammetric uses of Earth Observation data appear due to improvement of the geometric accuracy of the products. The use of a rigorous orientation geometric model for the georeferencing of satellite image data increases the geometric potential of the photogrammetric procedure relating the satellite image to the ground point coordinates. The application of the commonly used geometric models demands the use of several control points (GCPs); the exact number of GCPs depends on the unknown parameters that need to be defined. In this paper a new rigorous orientation geometric model for pushbroom linear array sensor based on collinearity equations is described. The model uses the interior orientation parameters and orbital data; in its full version the model has 24 unknown parameters. The advantages of this model in comparison with the normal use of collinearity equations are the existence of less unknown parameters and the use of three observation equations for each observed point on the image. The scale factor in collinearity equations is not eliminated but it is modelled as a second order surface polynomial.


INTRODUCTION 1.1 High Resolution Satellites
Most satellite agencies have been oriented to launch optical satellites that have the ability to acquire stereo pair images, usually in along-track, with high resolution; WouldView-1 & 2 and GeoEye-1 are commercial satellites with ground sample distance smaller than 0.5 m.Using array detectors which are arranged in a focal plane along a line (linear array sensors -pushbroom), with as long a width as possible, they can collect images with very high stability and give all parameters that are important to the photogrammetric processes, like focal length, principal point, orbital data, etc.Thus the use of a rigorous orientation model that allows increased accuracy in the photogrammetric process of these satellite images is very important.

Approximate and Rigorous Orientation Models
Many approximate models have already been developed for the geometric correction and the georeferencing of satellite images.Some of these are based on Rational Polynomial Coefficients (RPC), which are involved in metadata of images.The coefficients of polynomials have no physical substance.As it is well known, a Rational Function Model (RFM) is generally the ratio of two polynomials derived from the rigorous sensor model, ground control and/or sensor orientation information.These models do not describe the physical imaging process but they use a general transformation to describe the relationship between the image and the ground coordinates.For satellite sensors with a narrow field of view, like IKONOS, even simpler sensor models, such as the 3D affine model and the relief-corrected 2D affine transformation, can be used.Their validity and performance is expected to deteriorate as the area size and the rotation of the satellite during imaging (the so-called forward or reserve scanning mode), which introduces non-linearities, are increased.Also the 3D affine model cannot be applied when the height range is increased and when there is a poor Ground Control Point (GCP) distribution.Furthermore there are orientation models which are based on rigorous models, which use the physical imaging geometry.For pushbroom sensors in along-track direction, which means nearly parallel projection between lines and perspective projection in each line, the mathematical model is composed by the collinearity equations, including interior and exterior orientation parameters, combined with orbit determination-propagation models.The collinearity equation gives the geometry of a bundle of rays, describing the path of an optical ray from the object point on the ground to the sensor; moreover it gives the greater accuracy in order to relate two coordinate reference systems, for instance the ECEF (Earth CEntered Fixed) system and the camera's system; also, it provides for easier statistical analysis of the results.The mathematical solution depends on the sensor type; at any case GCPs are required and often many parameters are involved, some of them highly correlated; metadata files for satellite's orbit are required.The exterior orientation is modeled by piecewise polynomials, Lagrange polynomials or constant shifts.The interior orientation parameters are the focal length, the principal point displacements, the line curvature, the line rate, the lens distortion etc; some or even all of these parameters are not always involved in the mathematical model.

WorldView-2
The DigitalGlobe's WorldView-2 satellite was launched in October 8, 2009 and has the ability to collect panchromatic and multispectral images (in 8 bands) with pixel size of 0.46 m and 1.84 m at nadir respectively.WorldView-2 imagery products are delivered with a set of metadata files called Image Support Data (ISD).The proposed orientation model is applied only to Basic Imagery product, which includes Attitude, Ephemeris and Geometric calibration files.
The main technical characteristics of WorldView-2 sensor are: pushbroom camera with a focal length of 13246.139mm, image bands on Pan, Coastal Blue, Blue, Yellow, Green, Red, Red Edge, NIR1 & NIR2, 16.4 km swath and stereo capacity with Agile Telescope.

PROPOSED MODEL
The new rigorous model uses the collinearity equations combined with orbital and attitude data in order to relate the points in the ECEF object coordinate system to the corresponding points in the image coordinate system.The relationship between these two coordinate systems is based on: 1. three rotations, which are combinations of the Keplerian elements or the state vector, referred to the ECEF system, using the transformation parameters between the ECEF and ECI (Earth Centered Inertial) systems 2. three rotations (ω, φ, κ) for the additional undefined rotations of the satellite at the time of imaging and two off-nadir viewing angles (α, β) of the linear array sensor 3. the scale correction for each object point, between the two coordinate systems, expressed as a first or second order surface polynomial that depends on the ground relief.
This model can be implemented on pushbroom sensors (linear array sensors) with narrow Field of View (FOV).The scale correction, on account of narrow field of view and the significant difference between fly height of the satellite (770 km for WorldView-2) and ground elevation, can be modeled by a second order surface polynomial.The benefit of this model is that the individual point scale has not been eliminated and as a result one more observation equation added for each point; we have three equations for each measured image point instead of two.As a consequence we need fewer GCPs for determining the unknown's coefficients.So we have the possibility to figure out the exterior orientation for each linear array CCD in order to establish a relationship between the image points and the ground objects.
The observation equations of the rigorous orientation model yields: : image coordinate of the object point i S i (dt) : scale of the object i. ds(dx i ,dt) : scale correction of object i, expressed as a second order surface polynomial and is equal to: when dx i =(SensorColumns/2)-x i , with x i in pixel M q : the rotation matrix of pointing angles M b (dt) : the rotation matrix of the ω i (dt), φ i (dt), κ i (dt) angles around X, Y and Z axes, depending on dt for the i scan line M db (dt) : the rotation matrix of the dω(dt), dφ(dt), dκ (dt) angles, which are the corrections of the ω, φ, κ rotations: X gi , Y gi , Z gi : object coordinates of point i in ECEF system (WGS84) Xo(dt), Yo(dt), Zo(dt) : coordinates of the perspective center of the scan line, depending on dt dXo(dt), dYo(dt), dZo(dt): corrections of the perspective centers The model, in full version, has twenty four (24) unknown parameters (m 0÷5 , dω 0÷2 , dφ 0÷2 , dκ 0÷2 , dX 0÷2 , dY 0÷2 , dZ 0÷2 ) for the determination of the exterior orientation of each scene, which means that at least eight (8) GCPs must be visible on the whole image to figure out these parameters.Usually, the minimum number of the necessary GCPs is much smaller, since fewer parameters than those in the full version of the model are needed.The simulation and real data tests have proved that the use of only 10-13 parameters of the model depending on the ground relief and the distribution of the GCPs (see the following section) is sufficient; therefore in order to calculate the necessary parameters only 4 GCPs are needed.Consequently, the proposed rigorous model gets a significant advantage in comparison with the existing rigorous orientation models, especially for applications on large areas where the use of many scenes is needed.
Figure 1 shows the relation between ECEF (WGS84) and the camera's system according to the rigorous model.

Implementation of the Model
WorldView-2 images have been used for the implementation of the rigorous model.The selection of this sensor is made using as criteria the very small pixel size (0.5m) and the fact that it gives all necessary parameters for the implementation of a rigorous model like the collinearity equations.Relevant applications of the proposed model have also been made with other satellite sensors, e.g., using ASTER and QuickBird data, with very promising results in terms of the achived georeference accuracy and the elimination of the required GCPs.In the following an investigation of the existence of similar advantages in the most recent generation of satellites, with less than 0.5 m GSD, is made.
For the assessment, a bundle adjustment program has been performed using the MATLAB software environment that was developed in order to apply and check the proposed rigorous orientation model.Through the bundle adjustment the unknown parameters of the model were determined as well as the coordinates of the tie points.Figure 2 shows in flow chart the procedure of the bundle adjustment solution using the new rigorous model.
As mentioned the proposed model uses 24 parameters in the full version.For the presented research some versions of the model were evaluated, in terms of their ability to improve accuracies in both the object space and the image space.Also, it is necessary to examine the physical meaning, the significance and the contribution of each parameter in the exterior orientation procedure using statistical criteria.For example, the model might be more effective when of additional variables are included or perhaps when detection of one or more of its parameters is made.Adding a variable in the model always causes the sum of least squares to increase and the error sum of squares to decrease.We need to decide whether the increase in the regression sum of squares is large enough to justify the use of additional variables in the model.Furthermore, adding a statistically insignificant variable to the model may actually increase the error mean square, indicating that using this variable the model has poorer fit to the data.The hypotheses test would be useful in determining the potential value of each parameter of the model.After the performance of the bundle adjustment the calculation of the standard deviation, the covariance matrix, the correlation among parameters and the t-test parameters of the unknowns is made.If the standard deviation of the solution is better than 1 pixel, the solution is considered adequate.

WorldView-2 Data Set
The stereopair of WorldView-2 images used in this case study was acquired in May 4, 2010 in the Attica region, north-east of Athens, Greece (N 38.8 deg and E 23.7 deg), with product order id 10EUSI-1344-01_I002356_ST02-P005751.It is a suburban and semi-rural area (Figure 3) with semi-mountainous topography of 700 m elevation range (from 100 m to 800 m altitude).
The data set is a Basic Stereo Imagery product (Stereo 1B), panchromatic with pixel size of 0.55 m for the "Forward" image and 0.86 m for the "Reverse" image.The acquisition angles are -4 deg and -37.6 deg in track view for first and second taken images.They are supplied by the image support files .tif,.rpb,.imd,.geo,.eph,.att,.til,.ste.These metadata files are indispensable to extract information on interior orientation, camera geometry, acquisition time and exterior orientation.It is derived from the above data that the images used are rather a difficult case of WorldView-2 data, as they do not provide minimum possible GSD and the "reverse" image was taken with a great inclination, with obvious impact on the representation and the measurement accuracy of the objects (e.g., the buildings).
The determination of the unknown parameters was calculated having thirty three (33) GCPs and twelve (12) ICPs available.The proposed rigorous orientation model through a triangulation process was applied; thirty three (33) tie points (TPs) were measured.The ground coordinates of GCPs and ICPs were measured with GPS receivers, resulting to an accuracy of 10 cm in EGSA-87 (Greek national reference coordinate system).The image coordinates were measured manually.The distribution of GCPs and ICPs on the scene is illustrated in Figure 3.

Preparation of Auxiliary Data
For the application of the proposed rigorous model the auxiliary data that are supplied by ephemeris and attitude files (.eph, .att)into an appropriate type should be transformed.That means that the state vector of perspective center of camera must be in WGS-84, the attitude data that are given as quaternion elements must be transform in rotating angles in Euler type (ω, φ, κ) referring in X, Y, Z axes in WGS-84.The acquisition time of the first line scan for each image is then calculated.For the interpolation of the position and attitude, 3 rd order polynomials are used.The WorldView-2 satellite has the ability to estimate a state vector and attitude data every 0.02 sec by GPS and IMU instruments respectively, that are carried on the satellite.

Results
The first test is made using the full model, for each image of the stereopair, with 24 unknown exterior orientation parameters plus the focal length correction (df), in a simultaneous self-calibration adjustment.A significant test for each parameter is made using the t-test (student distribution) for confidence interval 95%.It is shown that only 14 unknown exterior orientation parameters are necessary in order to figure out the model without any degradation of the accuracy of the solution.Due to the semi-mountainous terrain it is proved that the use of a first order polynomial for the scale correction (a simple ruled surface) is sufficient.The fourteen significant parameters are: dX 0 , dX 1 , dX 2 , dY 0 , dY 1 , dY 2 , dZ 0 , dZ 1 , dZ 2 , dω 0 , m 5 , m 3 , m 2 , m 0 .Adding the df parameter for self-calibration, the total number of the unknowns becomes fifteen (15).So, the minimum number of GCPs for an independent adjustment of each image is only six (6).
Various solutions, using the full model or the statistically important parameters only, with a varying number of GSPs from the minimum required up to the maximum available, were made.Tables 1 and 2 show the residuals in image and ground space using characteristic cases of the bundle adjustment solution; the calculation of the X, Y, Z coordinates of the ICPs was made using intersections.

Assessment of the Results
Conclusions derived from the results of all tests can be summarized below:  Almost all solutions have achieved a sub-pixel accuracy in x and y.  Best results were achieved using the full model and the maximum number of GCPs.However, with only 15 parameters the accuracies in X and Y remains better than 0.5 m and in Z better than 1 m for the first image; the accuracy is worse mainly for X and Y for the "back" image, in which the acquisition angle is very large (about -37.6 deg) in track view and the pixel size is 0.86 m. show that the accuracy in elevations is rather independent from the number of GCPs, when we have more than 6 GCPs; the Y-coordinate is influenced more, although by using 10 GCPs the maximum accuracy that can be achieved is quarranted.The above results are derived from aerial-triangulation adjustments using the same 12 ICPs in all cases.

CONCLUTIONS
In this paper a new rigorous orientation model based on the collinearity equations is described, initiating a new method for the implementation of these equations.The proposed model can be applied for the orientation of linear array sensors using orbital data and sensor information.The benefit of this model is the usage of three observation equations, instead of two, for each measured image point, as the scale factor is not removed from the collinearity equations but is corrected by a second order surface polynomial as functions of dx i and dt i quantities.The full model has 24 parameters that all have physical meaning.In order to determine these parameters it is proved that at least six (6) GCPs are required, in comparison to the usual model of the collinearity equations, where at least nine (9) GCPs are required in order to calculate 18 unknown parameters.Application tests made using various simulation and real (spaceborn imagery) data have proved that usually fewer than 24 (of the full model) parameters are necessary for the successful application of the proposed model.
In this paper a stereopair of WorldView-2 satellite images was used with very satisfactory results.Applications have been done for almost all high resolution satellite sensors, e.g., ALOS, QuickBird, etc, with relevant results in terms of achieved (sub-pixel) accuracy and of the number of the statistically significant parameters of the model and the minimum required number of GCPs.In all cases 10-15 unknown parameters are sufficient for the achievement of the required accuracies; so, the number of necessary GCPs can be limited to six (6).
fl -t i , when t i is the time of the object point i or specific time of ephemeris data and t fl is the time of the first scan line f : focal length of the linear array imaging system.xo, yo : image coordinates of the principal point.y i

Figure 1 .
Figure 1.Relation between ECEF and the camera's system

Figure 4 .
Figure 4. RMS in X-coordinate of ICPs as function of GCP number

Table 1 .
Residuals of bundle adjustment with 25 unknown parameters

Table 2 .
Residuals of bundle adjustment with 15 unknown parameters