Direct georeferencing with correction of map projection distortions for active imaging

In aerial photogrammetry, the Cartesian coordinate system for the description of object space is commonly used. In contrast, many projects have to be processed in specific map projection and vertical datum. In that space, some geometric deformations exist. There are some compensative methods for active and passive sensors. In the case of active sensors, decomposition and the correction of observation vector for each ground point can be used. We obtain height, horizontal distance and horizontal angle in this process. All of these values should be corrected for precise georeferencing. The contribution deals with the derivation of the corrections and gets some theoretical values from the area of the Czech Republic. Esspetially in the case of high flying heights the corrections may gain values even in order of meters.


Introduction
The task of georeferencing in the field of aerial topographic survey is the determination of geometric relations between captured data and the real world [7].It includes two consecutive procedures:the determination of the exterior orientation parameters (EOP) and the restitution scene from EOP and observed data [9].EOP comprise of the spatial position of the sensor projection center and sensor attitude at the time of the observation.EOP for passive sensors can be gathered indirectly (by the measurement of image coordinates of ground control points) or directly from records of the on-the-board navigation system.On the contrary, active imaging is dependent (due to the sequential measurement principle andthe motion of the carrier vehicle) on direct methods.
Referring to the direct georeferencing EOP is typically measured by the GNSS and inertial measurement unit (IMU).GNSS provides absolute position with sufficient frequency (at least 1Hz) and IMU sensor attitude and acceleration.Final trajectory is deduced by combining GNSS and IMU observations [1].
The Aerial survey product is mostly provided in "national coordinates".It means the coordinate system realised by a combination of the national geodetic datum with a national map projection with the associated vertical system.The model for direct georeferencing is designed for cartesian space but national coordinates do not fulfill the condition and therefore cause various geometrical distortions.There are two ways to obtain accuratedata in national coordinates (See Fig. 1): restitution of the scene in cartesian (usually Earth-fixed) frame and transforming the scene to national coordinates or restitution of the scene in national coordinates with the corrected observation vector.Further text refers to active imaging, especially to georeferencing of LiDAR data.

Reference frames and EOP transformation
The EOP should be transformed into desired reference frame.Overview of reference frames purveys Tab 1. Detailed description of transformations can be found in [9] and [6].Brief summary contain Fig. 2.

Model for georeferencing of LiDAR data
According to [10], the coordinates of ground point X G can be (in cartesian frame) expressed as where X EO and R EO are sensor position vector and the rotation matrices formed by angular EOP, respectively, R scan and X range are scan angle matrix and range vector, respectively.The Geoinformatics FCE CTU 15(1), 2016 second term forms the observation vector for direct georeferencing X dg .If the georeferencing is carried out in national coordinates (p-frame), X EO can be given by exact formulas (usually provided by state mapping authority), contrariwise X dg is skewed due to datum scale change and p-frame geometry.We may assume that the correct position of ground point in p-frame X p trans G is given by georeferencing in e-frame and subsequently rigorous transformation to p-frame.Then we gain the correct observation vector in p-frame X pd g as Our task is to modify vector X p dg to X pd g .It will be to apply processes published in [10] as a practical approach.It involves some simplifications and refers to conformal map projection.

Correction of p-frame distortions
According to [10], the method consists of four subsequent steps (Fig. 3).

Datum scale distortion. The Cartesian e-frame and n-frame have a different scale (if they
have not used the same datum).Hence the length of the observation vector is different as well and it should be multiplied by scale factor m datum .
Decomposition of X p dg to height component Z dg (it is always negative), horizontal distance D and horizontal angle ϕ (Fig. 4).It will be called "spatial observations".IMU works based on Newtonian laws and its Z-axis is aligned with plumb line.If the vertical datum is not gravity-related, the gravimetric correction is required to be added in observations [5].Decomposition of X p dg is executed by simple equations

Application of map projection corrections
The Earth curvature correction Height and length component are deduce from geometry in Fig. 5 and Fig. 6. h ec can be expressed as G and F can be approximately regarded as locating at a same arc with the radius R + H G , then due to the small value of α we can simplify cos α≈1 where H S is the sensor projection centre.By combining ( 6), (7) we can obtain

Correction of length component
The first step is the conversion of the horizontal distance D to the geodetic distance S. S is always much shorter than the radius of reference sphere, therefore we can approximate the geodesic line by a circular arc.By using (7) we can form a relationship from the definition of the map projection scale factor we obtain the equation for calculation projected length where m stands for the map projection length factor.For this purpose it is sufficient to compute m in one point (preferably in the position of sensor projection centre).

Correction of horizontal angle
Since the normal-section-to-geodesic correction distinguishes negligible values, it should be applied to skew-to-normal and arc-to-chord corrections.Skew-to-normal corrections express the angle between the directions of the spatial straight line and its corresponding normal section.Regardless to used coordinate system it is given by [4] where ϕ G and h G are the geodetic latitude and ellipsoidal height of the ground point respectively, ρ m represents the curvature radius in the meridian plane and α is azimuth of observation vector.For practical calculations h G should be replace by h S + Z dg .
Arc-to-chord correction is angle between tangent of projected geodesic and its corresponding chord line.Computation is individual for each map projection.

Experiment
The asset of the above-mentioned procedure was proven in [10].Since the real LiDAR data is distorted by variety of random and systematic errors, the simulated data was used for the experiment.It was an applied method based on (2).It compared coordinates of ground points restituted in e-frame and transformed to p-frame (assume as the "correct" data) and the coordinates of ground points restituted in p-frame with application of correction.Most differences gain sub-millimeter values (Tab 2.)

Situation in the Czech Republic
The Earth curvature correction is area-independent thus we will discuss the other components in the most frequent national coordinates in the area of the Czech Republic: S-JTSK (Datum of Uniform Trigonometric Cadastral Network) / Bpv (Baltic Vertical Datum -After Adjustment) and UTM 33(34)N / ellipsoidal height (GRS80).Regardless to a projection used, we can compute the skew-to-normal correction as [3] ζ´= 0, 108 As the correction attains a value of 0,07 in extreme cases (Sněžka Mountain, azimuth 45°), we can disrespect it.

S-JTSK/Bpv
Change of the datum scale between WGS 84 to S-JTSK is -8.750 ppm, i.e if the observation vector is 1 km length it cause its shortening by 8.75 mm.This factor should be taken into account esspetialy in the case of high survey flights.
S-JTSK use double conformal conical Krovak projection (EPSG code 5514).Local scale is given by [2] where the first part marks constants α = 1.000597,N radius of the cross section and U spherical latitude.The projection has 2 undistorted parallels.The local scale causes distortion in the range from -10cm/km to 14cm/km.The difference between D and D in some selected locations shows Tab.3 Arc-to-chord correction is given (after removing high-order terms) as [8] where R means distance of the point from krovak s projection origin, D´i angle beetween X-axis and the point form projection origin and S i cartographic parallel.The undistorted parallel has been chosen as S 0 = 78 • 30´.For calculation parameters of the ground point, we can use the rough position obtained by the uncorrected vector X p dg .The magnitude of the correction depends on the length of the horizontal distance, orientation of X p dg and distance from undistorted cartographical parallels.It attains a size of few (in extreme case up to 10) arcseconds.

UTM/elipsoidal height (GRS80)
For practical applications, elipsoid GRS80 and elipsoid WGS 84 are identical, thus the m datum is the same.The local scale of each UTM strip is given by this simplified formula: where λ is longitude (measured from central meridian of the strip), ϕ is latitude and m 0 thescale factor of the central meridian.For the reduction of the scale distortion at the boundaries of the strips m 0 = 0, 9996 has been chosen.It follows the distortion -40cm/km at the central meridian and +17cm/km on the boundary of the strip.The length distortion of both projections shows Fig. 7.
After some simplification and removing high-order terms we can compute arc-to-cord correction as Geoinformatics FCE CTU 15(1), 2016

Conclusion
Mathematically, the only rigid way of direct georeferencing is to restitute the scene in an e-frame.However, after applying the above mentioned corrections, the residuals of the pframe deformations are negligible in comparison with the noise of GNSS/IMU measurements.
Although some formulas are simplified, remaining residuals are below 1 cm even by the relative flying height 8000m.The least impacted have angular corrections, and it can be ignored in the case of low relative flying height (up to 1000 m).
One of the main arguments for the choice of frames for georeferencing can be computational costs.The most demanding step comprises transformation of EOP to national coordinates (approximately four times higher computation time than the transformation of ground points).However, the GNSS/IMU record data with a frequency up to several hundred Hz and the stateof-the-art LiDAR systems emit pulses with frequency even 500 kHz [10].Therefore the number of ground points is much higher than the number of EOP observation and transformation of ground points cause consequently higher computation effort.
Georeferencing in national coordinates is strictly recommended to verify whether software take into account the p-frame corrections.

Figure 1 :
Figure 1: Accurate georeferencing in national coordinates for LiDAR.

Figure 3 :
Figure 3: Sequence of correction of observation vector in p-frame

Figure 5 :
Figure 5: Map projection distortion.G is the coresponding projected point of the G according to [10].

Figure 6 :
Figure 6: Detail for deriving the Earth curvature correction.O stands for centre of the reference sphere (rotated 90 degrees clockwise).

Table 1 :
Overview of the reference frames

Table 3 :
Difference between D and D in some selected locations.