EVALUATION AND ANALYSIS OF A PARAMETRIC APPROACH FOR SIMULTANEOUS SPACE RESECTION-INTERSECTION OF HIGH RESOLUTION SATELLITE IMAGES WITHOUT USING GROUND CONTROL POINTS

This paper initially proves mathematically that the solution of the relative orientation of the push broom stereo images leads to a set of ill-posed system of observation equations which means no relative orientation for this kind of imageries is possible. It is then demonstrated that the only way to bypass the inherent rank deficiency in the system of the relative orientation equations, is to incorporate fictitious data as known values into the equations. The proposed solutions given so far for the relative orientation of the push-broom images use a physical collinearity model and incorporate the so-called Virtual Ground Control points (VCPs) to bypass the ill-conditioning problem. The main objective of this paper is to replace the physical collinearity model with the 3D-affine transformation and use fictitious height information for the solution of the relative orientation of the push-broom images. The main advantage gained by this replacement model is its independence with respect to the sensor’s interior orientation parameters and also the simplicity of its functional model. Cartosat-1 P5 stereo images over a highly mountainous terrain are used to evaluate the potential of the 3D-affine for the solution of the aforementioned problem. The tests conducted in this research work indicate that the 3D-affine model with fictitious information can eliminate the ill-posed problem associated with the relative orientation of the push-broom images. However, the main disadvantage of the 3D-affine model is the fact that it imposes certain level of model deformation into the generated relative DEM. This paper also demonstrates the influence of the inclusion of the perspective-to-parallel transformation coefficient into the 3D-affine resection-intersection equations. The main contribution of this paper is, therefore, the solution of the relative orientation of the push-broom images with the 3D-affine model using fictitious information. The superimposed height deformation on the generated relative DEM is also analyzed.


INTRODUCTION
The so-called terrain independent Rational Polynomial Coefficients (RPCs) are generated by fitting the rational function models to the collinearity equations using precise attitude and positioning data generated by the sensors on board the spacecrafts.These supplied RPCs have opened a new era for the direct-geo-referencing technology.The RPC distortions, which are mainly confined to a simple shift for strip lengths of less than 50 km, can be compensated with a minimum of one ground control point (Fraser and Hanley, 2005).This means that the direct-geo-referencing has almost reached its final goal.
However, GCP measurement for RPC modification may undergo the following complications: (1) GCPs should be well defined points and identifiable on the ground.GCPs with this requirement are rarely available on satellite images that are acquired over nonurban areas, (2) GCPs must also be stationary, so that no displacement is occurred during the time lapse between the image acquisition and the field measurements, (3) even when the preceding conditions are fulfilled, the GCPs may not be accessible on the ground if they are located on highly steeped areas.In particular, for those projects that deal with monitoring dynamic phenomena such as land-slide, the above conditions are usually not fulfilled and consequently RPC modification may * Corresponding author fail.It is, therefore, important to find a reasonable alternative when the above mentioned complications arise.
One possible solution to tackle this problem is to generate a socalled relative DEM from the stereo satellite images without using GCPs.This relative DEM is then registered to the absolute DEM which is generated from the globally available elevation data.The so-called DEM-to-DEM registration (also called DEM matching) is then carried out by the seven parameters conformal transformation through which the relative DEM is connected to the absolute DEM using a curve matching approach.In this way, the complicated image-to-DEM registration is replaced by the process of DEM-to-DEM registration which is free from the complications and ambiguities associated with the image-to-DEM registration.
The above strategy does not require any real GCP and therefore the GCP complications mentioned above do not arise.However, the solution of the relative orientation of the push-broom images is impossible because the individual scan lines are disconnected and do not form a solid continuity and therefore, the system of equations for the relative orientation become ill-conditioned.One possible solution to deal with these situations is to incorporate a digital elevation model (DEM) generated from the globally available elevation data (e.g.SRTM) as a replacement model for the GCPs.However, the problem with these datasets is that the nodes of the DEM do not correspond to a definite feature on the image and therefore precise image-to-DEM registration may not be feasible.
The proposed solutions given so far for the relative orientation of the push-broom images use a physical collinearity model and incorporate the so-called Virtual Ground Control points (VCPs) to bypass the ill-conditioning problem, (Kim and Jeong, 2010).The main disadvantage of the physical models is the requirement for the sensor's interior orientation parameters which may not be available.Moreover, the algorithmic structure of the physical model is highly non-linear.This requires good approximations for the initial unknown values.In this paper, a parametric approach is implemented for the solution of the relative orientation of the push-broom images using VCPs.The aforementioned solution comprises three consecutive stages: -Relative orientation of the stereo satellite images (without using real GCPs), -Generating relative DEM by automatic or manual stereo-image matching, -Registration of the relative DEM to the absolute DEM.
The main objective of this paper is to implement and evaluate a parametric model for the solution of the first stage.Among the parametric models (e.g.DLT, 3D-affine, etc.), the 3D-affine is considered to be an optimum solution for handling the relative orientation of the high resolution push-broom images because this transformation conforms very well to the parallel nature of the high resolution satellite imaging geometry (Morgan et al., 2006).Although 3D-affine transformation has been widely in use for the geometric restitution of the push-broom images with the real GCPs, to our knowledge this transformation has not yet been tested and analysed for the relative orientation of the push-broom images with the virtual control points (VCPs).
In this paper, a 3D-affine transformation is implemented using VCPs for the simultaneous resection-intersection of the stereo Cartosat-1 P5 images acquired over a highly mountainous terrain.
A relative digital elevation model is then generated from the calculated height information.The main contribution of this paper is, therefore, the solution of the relative orientation of the push-broom images with the 3D-affine model using VCPs and also the analysis of the height deformation of the generated relative DEM by incorporating the perspective-to-parallel transformation coefficient in the resection-intersection model.In the section that follows, the impossibility of the relative orientation of the push-broom images is demonstrated mathematically.This is then followed by describing the process for generating VCPs.The solution of the relative orientation of the push-broom images with the 3D-affine transformation using VCPs and the evaluation of the precision of the generated relative DEM are then described.

GENERAL SCHEME OF THE PROPOSED SOLUTION FOR THE RELETIVE ORIENTATION OF THE PUSH-BROOM IMAGERIES
As mentioned in the introductory section of this paper, the relative orientation of the push-broom images is an ill-posed problem.The solution of the relative orientation for the images acquired by single projection cameras is feasible because in a frame camera, for any pixel located in the image frame, the x and y disparities caused by the physical rotations and translations of the camera are given by the following relations: Where ω,…, BZ represent the rotations and translations of the camera during the exposure time; and h denotes the height value on the ground.As can be seen from these equations, due to the rigidity of the frame cameras, y-disparities alone (Equation 1) may be employed to recover the full relative orientation parameters in a stereo pair without any recourse to the xdisparities.However, with the push-broom imageries, it can be shown that the sensor's rotations and translations generate x and y-disparities according to the following relations: Where the subscripts shift and drift denote respectively the constant and variable components of the translational elements for the push-broom sensors.As can be seen from Equation 2, ydisparities cannot be used to recover the full parameters of the relative orientation of a stereo pair.For the solution of the remaining relative orientation elements (i.e., BXshift, BXdrift and ) x-disparities must be used.However, the x-disparities, unlike ydisparities, are additionally influenced by the relief displacement.This means that the preceding parameters can be solved if and only if the relief displacement and consequently the corresponding height value on the ground is known for the tie points.Therefore, the only possible solution for this ill-posed problem is to include the object space information into the xdisparity equations.To avoid the GCP complications mentioned in the introductory section of this paper, the so called virtual control points (VCPs), as fictitious information, may be used instead of the real GCPs.
The general scheme in the form of a flow-diagram for the solution of the relative orientation (or simultaneous space resection-intersection) with fictitious information is depicted in Figure 1.Details of this diagram will be discussed in the sections that follow.
Figure1.Flow-diagram of the proposed method

THE PROCESS OF GENERATING VIRTUAL CONTROL POINTS
As expressed in the preceding section, the solution of the relative orientation leads to ill-posed set of observational equations for push-broom stereo images.The only way to tackle this problem without recourse to the GCPs, is to incorporate fictitious ground control information.However, this fictitious information should comply with the true terrain topography; otherwise, there will be deterioration in the precision of the generated DEM after the solution of the relative orientation.Nevertheless, if the fictitious height values are extracted from the grid nodes of the globally existing digital elevation data (e.g.SRTM), similarities between the generated relative DEM and the actual terrain topography will be preserved.However, the problem with these datasets, as mentioned in the introductory section of this paper, is that the nodes of the globally existing elevation data do not correspond to a definite feature on the image.
To overcome this problem, the greed nodes of the existing elevation datasets are transferred to the image space using approximate satellite ephemeris using the proposed method by (Kim et al., 2001).This transformation generates corresponding x and y pixel coordinates for any 3D grid node.The calculated pixel position and its corresponding 3D coordinates on the ground constitute a virtual control point (VCP).However, to incorporate these VCPs into the relative orientation equations, it is necessary to find a corresponding conjugate pixel on the second stereo pair.The position of this conjugate pixel can be determined by constructing a window around the pixel position of the left VCP and conducting one of the existing approaches of the image matching techniques.
It is important to note that these VCPs suffer from the misplacement errors.These misplacement errors, in their turn, cause datum displacement.Therefore, the datum plane of the generated DEM will be different from the initial datum and the absolute orientation for the generated DEM becomes necessary.
In the next section, the proposed approach for the solution of the relative orientation with the fictitious data is described.

PARAMETRIC SOLUTION OF THE RELATIVE ORIENATION WITH VCPS
As mentioned in the introductory section, the main objective of this paper is to implement the 3D-affine transformation for the solution of the relative orientation.The main advantages of this transformation are its simplicity and at the same time its compatibility with the imaging geometry of the high resolution satellite push-broom images.The 3D-affine transformation is given by the following relations: Where x, y, E, N, h denote the image and their corresponding object coordinates respectively; and A1 ... A8 are the 3D-affine parameters.However, the geometry in the pixel direction conforms to the projective geometry.The deviation in the pixel direction from the parallel geometry is directly proportional to the sensor's field of view.That is, the smaller the sensor's field of view, the smaller the deviation from the parallel geometry.Therefore, with the high resolution images the geometry in the pixel direction is very close to the parallel projection geometry.For a mountainous terrain, however, the aforementioned deviation is amplified.In this case, the projective geometry in the pixel direction can be made parallel, with the so-called perspective-to-parallel transformation.Thus, to reduce the model deformation associated with the relative orientation of the 3Daffine method, a perspective to parallel transformation coefficient, which expresses the viewing angle of the pixels as a function of the relief variations on the ground, is incorporated in the adjustment model as indicated by the following relation for the along-track stereo images (Morgan et al., 2006): To solve simultaneous space resection-intersection using these equations the linearization is required.The linearized form of equation ( 4) are given by the following relations:  =  1 +  2 + ℎ 3 +  4 +  1  +  2  +  3 ℎ (5)  =  5 +  6 + ℎ 7 +  8 +  5  +  6  + ( 7 +  tan )ℎ The above linearized equations are written for the selected VCPs and the well distributed tie points over the overlap area of the stereo pair and are solved iteratively in a simultaneous resectionintersection adjustment.With n tie points and m VCPs, there will be 4(n+m) equations in 16+3n unknowns.Therefore, to avoid illconditioning problem, tie and VCP points must satisfy the condition 4(n+m) ≥ 16+3n.The design matrix, constructed for 8 VCPs and 9 tie points, assuming unit weight matrices, is presented schematically in Figure 2. The sub-matrices depicted in the left side of the design matrix (Figure 2) are the partial derivatives with respect to the parameters of the 3D-affine transformation and the right side submatrices are the partial derivative with respect to the 3D object coordinates of the tie points.The advantage of the transformation presented in Figure 2 is its simplicity and its independence with regard to the sensor's interior orientation parameters.
Solving the linearized equations for VCPs and tie points, the 3D affine parameters for each stereo pair are computed.Using these parameters and dense matching for the entire overlap area, the 3D-afine space intersection will be solved to create the 3D ground data.A digital elevation model is then generated based on these data.As mentioned before, these 3D ground data are derived from the calculated elements of the relative orientation described in the preceding paragraphs which means that the generated DEM is relatively oriented.Therefore, this process should be followed by performing the absolute orientation with the seven parameters similarity transformation.This similarity transformation will then transfer the relative DEM into the absolute DEM provided that model deformation has not been occurred.If model deformation exists then the seven parameters conformal transformation may be replaced by a twelveparameter affine transformation during the process of DEM-to-DEM matching.This issue is a topic for future investigation.

EXTPERIMENTAL EVALUATION OF THE PROPOSED STRATEGY FOR THE SOLUTION OF THE RELATIVE ORIENTATION
Cartosat-1 P5 stereo along track images acquired over a highly mountainous terrain (Roodehen, Iran) are used to assess the proposed strategy for the solution of the relative orientation.The stereo images in reduced resolution are presented in Figure -3.The specification of the sensor is given in Table 1.The pixel coordinates of the VCPs, as discussed before, is calculated by an inverse transformation from the object space onto the forward Cartosat-1 image using the ephemeris provided by the metadata.The corresponding pixel position on the forward image is determined using manual stereo-matching.With these VCPs the relative orientation of the stereo pair is performed.The solution of the system of equations converged with 4 iterations.Therefore, the test results clearly indicated that the illconditioning problem mentioned before, can be handled well with the VCPs.Automatic image matching for the entire overlap area of the stereo pair is then carried out.Using space intersection the 3D object coordinates and subsequently the relative DEM is derived.
To assess the precision of the generated relative DEM, supplied RPCs from the Cartosat-1 metadata are used to generate 3Dobject coordinates.These RPC-based generated object points have high level of precision (Fraser & Hanley, 2005) and therefore are used to assess the precision of the generated relative DEM.This is achieved by calculating the object coordinates for all of the conjugate pixels using RPC-intersection.These 3D object points are then compared with the object coordinates of the same conjugate points derived by 3D-affine intersection.The scatter plot of the height residuals versus E-coordinate is give in Figure 6(a).The RMSE of the residual height errors, as depicted in Table 2, highly exceeds the expected values.This is due to the fact that the pixel positions of the VCPs are in the featureless areas and as a result, manually determined conjugate points suffered from the misplacement error.To improve the relative orientation precision, least squares image matching was conducted to find the conjugate points for the VCPs with higher level of precision.

RMSE of height errors (Meters)
Manual stereo-matching 16.85 Least squares matching 10.04 Table 2. Effect of VCPs position improvement using LSM The RMSE of the residual errors for the 3D components of the tie points in the object space after performing least squares image matching for the VCPs are given in Table 2.The corresponding scatter plot of the residuals for the tie points are presented in Figure 6(b).The results indicate improved precision as compared with the manual approach for the VCPs conjugate point determination.However, as the scatter plots indicate, the random components of the residual errors still are larger than what is expected from the 3D-affine transformation.The sources of these amplified random components are partly related to the errors in stereo matching for the VCPs and partly due to the imprecision of the original elevation data and also to some extent due to the VCPs misplacement errors in the image space.
In addition to these random components, the scatter plots of the residual height errors, presented in Figure 6, also indicate a systematic trend.This systematic trend shows a non-linear behaviour as presented in Figure 6 To compensate such a model deformation, a higher order transformation for the absolute orientation is required.
For the assessment of the nature and the contribution of each of the aforementioned factors which influence the precision of the relative orientation and consequently the precision of the generated relative DEM, further investigation on different types of terrain topography is being conducted in our institute and will be reported in future.

DISCUSSION AND CONCLUSION
In this paper, a 3D-affine transformation for the solution of the relative orientation or space resection-intersection of the pushbroom stereo pairs with virtual control points is implemented and analyzed.The initial tests conducted over highly mountainous terrain topography indicate that the inclusion of the VCPs can resolve the ill-posing problem associated with the relative orientation of the push-broom images.The proposed parametric approach for the solution of the relative orientation in some respects has superiority over the conventional physical approach because it does not require the interior orientation parameters and also due to the simplicity and linearity of the functional model wider range of initial approximations for the unknowns is required as compared with the physical models.Nevertheless, due to the intervening error sources, the precision figures obtained over the test area are not satisfactory.The main error sources can be summarized as follows: -Errors due to the misplacements associated with the image coordinates of the VCPs, -Height errors associated with the grid nodes of the incorporated globally available DEM, -Errors associated with the 3D-affine transformation, in particular the errors related to the deviation of the perspective projection from the parallel geometry in the scan line direction.
The first source of the error, introduces amplified random components and systematic error depending on the nature of the VCPs misplacements.The Second and third error sources generate systematic errors leading to the 3D model deformation for the generated relative DEM.The amount of the introduced model deformation also to some extent depends on the terrain topography undulation.That is, in highly oscillating terrain surface, the VCPs misplacements are more likely to engender higher level of the model deformation.In a terrain with a smoother level of the height undulation, on the contrary, the generated relative DEM is less susceptible to systematic deformation errors.In the former case, the remaining model deformation may be corrected during the process of relative to absolute DEM registration by replacing the seven parameters conformal transformation with the higher order transformations.Further analysis of the error sources and the strategies to deal with the model deformations should be conducted in future investigations.Moreover, the inclusion of the weight matrix in the simultaneous resection-intersection equations may be beneficial as regards the reduction of the degree of the influence of the image coordinates misplacement of the VCPs in the final results.This also requires further investigations.

Figure 2 .
Figure 2. Structure of the design matrix

Figure 3 .
Figure 3. CARTOSAT-1 stereo pair over Roodehen, IranVirtual control points for the solution of the relative orientation are generated using SRTM data for the area covered by the stereo images.The shaded relief of the SRTM data is presented in Figure4.

Figure 4 .
Figure 4.A view of the shaded relief of the Roodehen, generated from the SRTM data The proposed methodology, depicted by the flow-diagram shown in Figure 1, is implemented in MATLAB software.As the flowdiagram shows, the workflow begins with selection of VCPs on the SRTM data.As mentioned before, the minimum number of the VCPs for the solution of the relative orientation is 8, however, to increase the redundancy, 12 VCPs are considered for the solution of the system of equations.Consequently, 12 well

Figure 5 .
Figure 5. Distribution of VCPs over SRTM data Figure 6.(a): Scatter plot of the height residuals when VCPs conjugate point determination performed manually, (b): when VCPs conjugate points derived by least square image matching.
(a), whereas the height residuals given in Figure6(b) indicate a rather linear systematic trend.This linear trend is possibly due to the fact that the misplaced positions of the VCPs in the image space has introduced a shift and rotation into the calculated parameters of the relative orientation causing a corresponding shift and rotation in the datum plane.As mentioned earlier, this datum shift and rotations can be compensated by performing the seven parameters conformal transformation by which the relative DEM is transferred to the plane of the original absolute DEM.However, the non-linear behaviour of the residual trends, Figure6(a), is more indicative of a model deformation introduced by the non-linear VCPs misplacements in the image space.

Figure 7 (
Figure 7 (left) shows a view of the reduced resolution shaded relief of the generated relative DEM.The shaded relief of the object points derived by the RPC-intersection is given in Figure 7 (right).