GEODETIC REFERENCE SYSTEM TRANSFORMATIONS OF 3D ARCHIVAL GEOSPATIAL DATA USING A SINGLE SSC TERRASAR-X IMAGE

Many older maps were created using reference coordinate systems which are no longer available, either because no information to a datum was taken in the first place or the reference system is forgotten. In other cases the relationship between the map’s coordinate system is not known with precision, meaning that its absolute error is much larger than its relative error. In this paper the georeferencing of medium-scale maps is computed using a single TerraSAR-X image. A single TerraSAR-X image has high geolocation accuracy but it has no 3D information. The map, however, provides the missing 3D information, and thus it is possible to compute the georeferencing of the map using the TerraSAR-X geolocation information, assembling the information of both sources to produce 3D points in the reference system of the TerraSAR-X image. Two methods based on this concept are proposed. The methods are tested with real world examples and the results are promising for further research.


INTRODUCTION
The integrated use of three dimensional Archival Geospatial Data (AGD) requires an established relationship between the data and a Global Geodetic Reference System (GGRS).GGRS 1) ensures the geometric interoperability of data (existing and new), 2) it facilitates inter-temporal worldwide seamless mapping and 3) it is imposed in any use of the GNSS technology either directly (for instance in situ measurements) or indirectly (for instance air-born/space-born photogrammetry).Normally, the relationship between AGD and GGRS is established indirectly through the relationship between AGD and a regional geodetic reference system.Then the national definition of the regional geodetic reference system is used in order to convert AGD to GGRS.However in some cases, especially when it comes to older AGD, the conversion is not always a straightforward process: the definition of the regional geodetic reference system may be unknown, lost, ambiguous or it may not be readily available at the time of processing.In these cases the recovery of the relationship between AGD and GGRS is a laborious process as it requires in situ measurements.Alternatively, other 3dimensional data (mostly heterogeneous multitemporal and multimodal) with established relationship with GGRS may be used for this purpose, depending on the project standards (accuracy, time, budget, type of available data etc).It is not rare that none of these two alternative options is available or possible.This paper proposes a third option, the recovery of the relationship between AGD and GGRS by inverting the resection in space between a single, 2dimensional, SSC TerraSAR-X image and the AGD.The process takes advantage of innovative characteristics of modern satellite SAR sensors: 1) the high resolution imaging technology, 2) the accurate and ground-independent georeferencing capabilities and 3) the rich metadata that accompanies each image (Moreira et al, 2005), (Fritz et al, 2007), (Yoon et al, 2009), (Vassilaki et al, 2011).In essence control points on the image are "projected back" to the object space using height information of the AGD.
The resulting 3D global coordinates of the control points, together with the map local coordinates, can then be used to recover the relationship between AGD and GGRS.Two different methods are proposed in this paper: • The indirect method uses the GGRS and the AGD local coordinates of the control points to compute the parameters of the relationship, employing the Least Square Adjustment (LSA) .While for many relationships the LSA is linear, the "back projection" is not linear and requires a reliable first approximation.
• The direct method computes the parameters of the relationship directly from the AGD local coordinates and the Terra-SAR-X image coordinates of the control points, employing the LSA.In this case the LSA is highly non-linear and needs a reliable first approximation, which is taken from the results of the indirect method.
Real world data over two different sites is used in order to apply and evaluate the proposed process.Results are encouraging for further research.

INVERTING THE SAR PROJECTION
The most important part of the georeferencing using a single TerraSAR-X image is to produce 3D points in the object space (GGRS), using image coordinates and elevations from the map, or in other words, to invert the SAR projection.The following description of the TerraSAR-X projection is taken from (Vassilaki et al, 2011), (Vassilaki and Stamos, 2014) where more details can be found.

The SAR projection
The SAR projection transformation projects the 3D object coordinates of a point P (XP , YP , ZP ) to its 2D projection point p(xp, yp) on the SAR image.On the image, yp is the translated and scaled time stamp tp when the sensor measured point P, and xp is the translated and scaled distance RP between point P and Figure 1.SAR imaging geometry.
the sensor S(XS(tp), YS(tp), ZS(tp)) at the time of the measurement tP (Figure 1): The four translation and scale parameters R0, ∆R, t0, ∆tAz are taken (or computed) from the rich metadata that accompanies the TerraSAR-X image (Fritz et al, 2007), (Vassilaki et al, 2011).The distance RP is computed as the distance of point P and the location of the sensor at the time of the measurement.The distance is given by: (2) The timestamp tp is given by the solution of the following equation: where the coefficients ki of the polynomial are functions (Vassilaki et al, 2011) of the coefficients ai, bi, ci of the orbit of the sensor: The coefficients ai, bi, ci are computed using the locations and the velocities of the sensor at two time stamps which enclose the duration of the acquisition of the TerraSAR-X image and are provided by the rich metadata which accompanies the SAR image (Vassilaki et al, 2011).

Additional equation
Clearly the projection transformation (equation 1) can not be inverted unless an additional equation is used.This equation can be provided by the orthometric elevation H which is taken from the map.
where h is the geometric elevation and ∆h is the geoid undulation at the centre of the SAR image.Undulation varies very little compared to the accuracy of the elevation, and thus it can be considered constant without loss of accuracy.The computation of geometric elevation h is more complicated and depends on the geocentric coordinates of point P and the parameters of the reference ellipsoid: where a is the semi-major axis and b is the semi-minor axis of the reference ellipsoid.

Computation of the inversion
Equations 1 and 5 form a system of 3 equations with 3 unknowns XP , YP , ZP : where RP , tP , h are highly non-linear functions of XP , YP , ZP .The system can be solved by the non-linear Newton-Raphson iterative numerical method (Press et al, 1992).The method needs initial values for the geocentric coordinates XP , YP , ZP of point P .Luckily the XML metadata file which accompanies the TerraSAR-X image (Fritz et al, 2007) contains the geodetic coordinates λ, φ of the centroid of the image and the average geometric elevation h of the image.These can be found in the following section of the XML file: • productInfo/sceneInfo/sceneCenterCoord/lat • productInfo/sceneInfo/sceneCenterCoord/lon • productInfo/sceneInfo/sceneAverageHeight and they can be converted to geocentric coordinates (Redfearn, 1948).Extensive numerical experimentation showed that the coordinates of centroid are sufficient to trigger convergence of the method for any point of the SAR image.
The method also requires the derivatives of the equations with respect the geocentric coordinates XP , YP , ZP , which are rather involved and are shown below.The derivatives of the range (xp) image coordinate are: The derivatives of the azimuth (yp) image coordinate are: where: and where: A closer examination of equations 6, show that the computation of the geometric elevation h is in reality numerical.Thus the derivatives of the geometric elevation h are also computed numerically (Press et al, 1992).
The computation of the geoid undulation is done using the Earth Gravitational Model 2008 (EGM2008) (Pavlis et al, 2008), which is considered to be among the most accurate global geoid models according to the evaluation of the International Centre for Global Earth Models (ICGEM).

COMPUTATION OF THE GEODETIC TRANSFORMATION
The relationship between the GGRS and AGD local coordinates is the similar transformation, as the relationship is comprised of a translation and a rotation.A scale factor is also considered in order to take into into account any contractions or expansions of the paper, as most AGD is in paper form.The similar transformation is 2D as the SAR image lacks elevation information, and in fact the elevation is taken by the AGD in order to compute the GGRS coordinates of the control points.The 2D similarity transformation implies that the geocentric coordinates of the control points must be transformed to 3D coordinates of a geodetic projection such as the Transverse Mercator.
Two methods for the computation of the parameters of the similarity transformation are presented, the indirect and the direct method.It is trivial to extend the indirect method to compute the parameters of other relationships, while this is more involved for the direct method.

Indirect method
The similarity transformation considers two translations Xo, Yo, a rotation φ and a scale µ between the GGRS X, Y and AGD local x, y coordinates: The transformation is non-linear but can be turned into a linear one by setting a = µ cos φ and b = µ sin φ: In the indirect method, equation 8 is applied to the GGRS and AGD local coordinates of all control points.The GGRS coordinates of the control points are computed by inverting the SAR projection as described in Section 2. The parameters Xo, Yo, a, b are computed employing the linear LSA.

Direct method
In the direct the method, the similarity transformation and the TerraSAR-X projection are applied simultaneously to the AGD local coordinates of the control and yield the TerraSAR-X image coordinates.Specifically equation 8 is used to transform AGD local to GGRS coordinates, and then equation 1 in conjuction with equations 2-5 are used to compute the TerraSAR-X image coordinates from the AGD local coordinates.Combined, all the equations can be written as: The direct method is superior to the 2-step indirect method, as potential bias present in the similarity parameters will be propagated to the computation of the TerraSAR-X image coordinates and thus the LSA will tend to correct it.
Figure 2. Distribution of the GCPs on the TerraSAR-X image in Attiki.

APPLICATION
The method was tested with real world data over two different test sites in Greece: i) the Attiki and ii) the Thassos test site.

Test site in Attiki
The test site in Attiki is located in the greater north-eastern region of Athens.The area has steep mountainous terrain and it is generally covered by sparse vegetation.Kalamos and Markopoulo Oropou, are two communities localised in the area.The TerraSAR-X image spans almost entirely 2 maps and a small part of other 2 maps.Thus results using 1, 2 or 4 maps were computed.All maps contain a well defined grid which can be used to correct maps from defects due to distortions of the paper.
Thus results with and without corrections were also obtained.
Figure 4.A trigonometric point on the map which is used as CP.

Control and Check Points.
In order to apply the proposed method 29 Ground Control Points (GCPs) were identified both on the maps and on the SAR image (Figure 2).An effort was made to distribute the GCPs so that they cover the whole scene of the SAR image and the maps.This proved to be difficult as the map is much older than the SAR image.Almost all GCPs are located on crossroads which makes their position somewhat ambiguous.Some characteristic points of construction, such as building corners and fences, were also identified, but almost all proved to be unreliable, as most of the constructions had been rebuilt (Figure 3).
The maps also contained the locations of trigonometric points, the coordinates of which are known independently in the Hellenic Geodetic Reference System (HGRS87  specifically to the map which is covered almost entirely by the TerraSAR-X image.Then the maps were corrected using the grid points available in the maps.The grid points of the maps are aligned and thus, once corrected, 2 or more maps can be combined into a single bigger map.The methods were applied to combined maps comprised of 1, 2, or 4 maps. The geocentric coordinates of the GCPs, which are computed by both methods, are transformed to HGRS87 in order to match the reference system of the trigonometric points.The Root Mean Square Error (RMSE) of the transformation computed with CPs is given in Table 1 for all test cases.The LSA fit error computed with GCPs is shown in Table 2.The detailed errors of the trigonometric points used as CPs for a map comprised of 2 corrected maps are shown in Tables 3 and 4 for the indirect and the direct method respectively.
Overview of the tables show that while the indirect and direct methods produce different results, the differences between them are negligible.Thus one might prefer the indirect method as simpler.
The RMSE for the uncorrected map is larger than that of the corrected map, perhaps due to distortions introduced when it was photocopied from the original blue prints.The RMSE for 1 corrected map is 3.6 m which is good if one considers the huge temporal difference between the maps and the TerraSAR-X image.
The RMSE tends to increase slightly as the number of the maps increases.This is possibly because the similarity transformation may not accurately model the conversion from the system of the map (Azimuthal Projection) to the system of HGRS87 (Transverse Mercator projection) in extended Finally it is noted that RMSEs shown in Table 1 correspond to 4-5 pixels in the SAR image.

Test site in Thassos
The test site is in Thassos island of northern Greece, in the vicinity of the town of Thassos.The area has steep mountainous terrain and it is generally covered by dense vegetation.The town of Thassos is localised in the area.

Control and Check
Points.In order to apply the proposed method 13 GCPs were identified both on the maps and on the SAR image (Figures 5, 6).An effort was made to distribute the GCPs so that they cover the whole scene of the SAR image and the maps.This proved to be impossible as the map is much older than the SAR image.Again almost all GCPs are located on crossroads which makes their position somewhat ambiguous.
The maps also contained the locations of trigonometric points, the coordinates of which are known independently in HGRS87.
The trigonometric points were used as Check Points in order to test the accuracy of the methods.

4.2.4
Results.The georeferencing of the maps was computed using both the proposed indirect and direct methods.The methods were at first applied to one uncorrected scanned map.Then the maps were corrected using the grid points available in the maps, and the methods were applied to one corrected map and a combined map comprised of 2 maps.
The geocentric coordinates of the GCPs which are computed by both methods are transformed to HGRS87 in order to match the reference system of the trigonometric points.The RMSE of the transformation is given in Table 5 for all test cases.The LSA fit error computed with GCPs is shown in Table 6.The detailed errors of the trigonometric points used as CPs for a map comprised of 2 corrected maps are shown in Tables 7 and 8  It must be noted that RMSEs shown in Table 5 correspond to almost 2 pixels in the SAR image and in Attiki test site the RM-SEs correspond to 4-5 pixels.This means that, given the accurate and ground-independent georeferencing capabilities of the TerraSAR-X sensor, the RMSE is dominated by other sources of error.are a) absolute error of the map which is 2.5 m horizontally and 4 m vertically b) the huge temporal difference between the maps and the SAR image, which makes the placement of control points more or less unreliable, and c) the fuzzy and speckled nature of the SAR image, which makes the identification of control points ambiguous.Use of linear features instead of points, would alleviate the errors due to b) and c) and thus it would probably increase the accuracy significantly.

CONCLUSIONS
Two new methods for the georeferencing of AGD (maps) using a single TerraSAR-X image were presented.The methods exploit the accurate and ground-independent georeferencing capabilities of TerraSAR-X images and the elevation information of the map in order to invert the SAR projection transformation and produce 3D points in the object space, which are used for the computation of the georeferencing.Alternatively, the georeferencing is computed directly in a single step.The methods are cost effective as they are able to do the georeferencing without the obligation to use a pair of images and the necessary lengthy process to create a DTM out of them.The methods were tested with real world data in two different test sites and give satisfactory results.For further research the use of GCLFs (Ground Control Linear Features) instead or in combination with GCPs is proposed as GCLFs show good performance in multimodal and multitemporal image analysis.The use of the forthcoming TanDEM-X global DEM may also enhance further the process.
Yo, a, b; x, y, H)    yp = f2(Xo, Yo, a, b; x, y, H)(9)where xp, yp are the TerraSAR-X image coordinates, x, y are the AGD local coordinates, H is the orthometric elevation obtained from AGD, Xo, Yo, a, b are the parameters of the similarity transformation, and f1, f2 are highly non-linear functions which encompass equations 8 and 1-5.The system of equations 9 can be solved for the unknown parameters Xo, Yo, a, b by the non-linear LSA.The non-linear LSA requires a reliable first approximation of the similarity parameters, which can be obtained by employing the indirect method, before employing the direct method.The derivatives of functions f1, f2 of equations 9 with respect to x, y, H are also required by the non-linear LSA.Under the chain rule the derivatives are comprised of the derivatives of equations 8 and 5 which are trivial given that the geoid undulation is practically constant, and the derivatives of equations 1 which are rather involved and are given in Section 2.3 4.1.1TerraSAR-X data.The TerraSAR-X image used in the test is a basic image product (Figure2), captured on February 2009 with the High Resolution SpotLight (HS) acquisition mode and it is of type Single Look Slant Range Complex (SSC).The incidence angle is 53 • and the polarisation is HH.The projected spacing values for range and azimuth are 0.45 m and 0.87 m, respectively.

Figure 3 .
Figure 3.Some homologous GCPs on the map (top) and on the TerraSAR-X image (bottom).
Figure 5. Distribution of the GCPs on the TerraSAR-X image in Thassos.
-X data.The TerraSAR-X image used in the test is a basic image product (Figure5), captured on March 2011 with the High Resolution SpotLight (HS) acquisition mode and it is of type Single Look Slant Range Complex (SSC).The incidence angle is 41 • and the polarisation is HH.The projected spacing values for range and azimuth are 0.91 m and 1.81 m respectively, which are double of those in Attiki test site.4.2.2 Archival data.The archival data (AGD) which was used in the study belong to the same series as in the Attiki test site (Figure6).The TerraSAR-X image spans 2 maps.Thus, results using 1 or 2 maps were computed.

Figure 6 .
Figure 6.Distribution of the GCPs on the maps in Thassos.

Table 1
). HGRS87 uses a Transverse Mercator map projection and is based on a shifted version of the GRS80 ellipsoid.The trigonometric points were used as Check Points (CPs) in order to test the accuracy of the method.
. The RMSE of all test cases in Attiki test site, computed using CPs (trigonometric points).

Table 2 .
The LSA fit error of all test cases in Attiki test site.

Table 3 .
The errors of CPs for the indirect method applied to 2 corrected maps in Attiki test site.

Table 4
. The errors of CPs for the direct method applied to 2 corrected maps in Attiki test site.

Table 5 .
The RMSE of all test cases in Thassos test site, computed using CPs (trigonometric points).

Table 6 .
The LSA fit error of all test cases in Thassos test site.

Table 7 .
The errors of CPs for the indirect method applied to 2 corrected maps in Thassos test site.

Table 8 .
The errors of CPs the direct method applied to 2 corrected maps in Thassos test site.andthe direct method respectively.Overview of the tables show again that the indirect and direct methods produce hardly different results.The RMSE of the corrected map practically matches that of the corrected map, which means low distortion of the paper and/or that other error dominates the RMSE.The lower RMSE for the 2 maps is attributed to the bigger number of CPs.The RMSEs in this test site (Thassos) are generally larger than the ones in test site Attiki, probably to the lower resolution of the TerraSAR-X image.