Determination of geographic position of object by applying 3D polar observations/Objekto geografines padeties nustatymas taikant erdvinius polinius matavimus.

1. Introduction

The unambiguous position of objects on the surface of the earth or in the space around it can be determined in the selected solid geodetic coordinate system.

The position of an object in a space of local rectilinear Cartesian coordinates is simply determined by 3D polar observations. A simple geometric observation model and computations can be used in such a space, but application of such polar observations in a space of the global geographical coordinate system is complicated, because we come into collision with the curvilinear coordinate system related to the surface of an ellipsoid. We get a complex geometric observation model in such case.

There are known methods to determine the position of an object by 3D polar observations when the coordinates of the measurement instrument are well-known. There are no methods or algorithms to determine the coordinates of measurement instruments by polar observations to a point with well-known coordinates, however (e.g., to determine the position of a measurement instrument in an aircraft by measurement to a point with well-known coordinates).

The objective of this work is to substantiate a technological scheme and computational algorithm which simplifies the solution of the problem for determining geographical position of an object with the extension of this algorithm for the solution of the inverse problem when the position of the measurement instrument may be determined by measurement to a point with well-known coordinates.

Simplification can be achieved by applying an auxiliary local 3D Cartesian horizontal coordinate system. In that way, the complicated measurement model linked to the variable curvature of earth's ellipsoid surface in 3D used for the determination of object geographic placement using the polar method is avoided.

2. Geometric model of space-polar observations

The geographic position of an object is described by geodetic coordinates: B--geodetic latitude, L--geodetic longitude, and [H.sup.(g)]--geodetic height. The LKS94 national geodetic coordinate system that is presently used in Lithuania is the realisation of the ETR89 European geodetic coordinate system in Lithuania (Skeivalas 2008; Zakarevicius 1994; Zakarevicius 2003). The WGS84 international geodetic coordinate system is used for international geodetic and cartographic work under NATO requirements (Boucher, Altamimi 2001; Logan et al. 2003). According to the requirements of International Civil Aviation Organization (ICAO), the WGS 84 coordinate system is used in aviation (the ETRS 89 coordinate system may be used in Europe). Therefore 3D polar observations for the determination of the geographic position of objects should be relative to these geodetic coordinate systems.

For determination of object spatial position by the polar method, three quantities must be observed: distance from the origin of polar coordinate system to the object to be determined and two angular quantities. These quantities cannot be observed directly in the geodetic coordinate system. Therefore an auxiliary coordinate system is needed for direct polar observations. This auxiliary coordinate system may be used for easy and unambiguous relationship with geodetic coordinate system. Such an auxiliary coordinate system may be a 3D Cartesian horizontal coordinate system with its origin at the position of the measurement instrument, Geodetic or 3D Cartesian geocentric equatorial coordinates of point should be known for relating observations to the geodetic coordinate system. Plane xoy of the horizontal coordinate system passes through observation equipment and is perpendicular at this point to normal of the ellipsoid describing the geodetic coordinate system (Fig). The initial point [P.sub.o] may be arranged at a fixed place or in a moving object (e.g., aircraft). In the case of point [P.sub.o] in a moving object, its coordinates must to be determined at the moment of measurement. This may be done, for example, by using the global positioning system (GPS).

Coordinate plane xoz is the meridian plane at the origin of the coordinate system. Coordinate plane yoz is a prime vertical plane at the origin of the coordinate system. The positive direction of the x axis of the horizontal coordinate system points to the North Pole, y points to the east, z points to the zenith, i.e. z axis corresponds at the point of origin with the normal to ellipsoid.

To determine the 3D position of object [P.sub.i], both distance [s.sub.oi] and the horizontal angle are measured from the meridian plane passing through the coordinate origin to the normal plane passing through origin and the object to be determined, i.e. the direction from the origin to object [A.sub.oi] and the zenith distance (angle) to the object to be determined [z.sub.oi] (Fig).

[FIGURE OMITTED]

The 3D Cartesian horizontal coordinates of object [P.sub.i] will then be

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

The position of a 3D point in the local horizontal coordinate system with the origin coinciding with instrument location is determined in this way.

This geometric model of polar observations is suitable for the geodetic horizontal coordinate system (when the z axis coincides with normal direction to ellipsoid). It can be also applied for coordinate system of astronomic horizontal when z axis coincides with vertical (gravity) direction.

Instruments of observation during measurement are oriented not along the normal but along the vertical (gravity direction). Direct measurements are therefore done in the astronomical coordinate system but not in the geodetic coordinate system. The 3D position of the objects should however be defined in a geodetic coordinate system. Corrections to the measured values of azimuth and zenith distance should therefore be computed. It is necessary to know an angle between Normal and Vertical (vertical deflection) at the instrument position for correction computation. Vertical deflection is described by two components: [xi]--value of vertical deflection in the meridian plane and [eta]--vertical deflection value in the first Vertical plane. If vertical deflection values at the origin of the horizontal coordinates are [[xi].sub.o] and [[eta].sub.o], the corrections for the observed azimuth and zenith distance will be (Jakucionis et al. 1999; Moritz 1990; Zakarevicius 1998):

[[DELTA].sub.A] = [[eta].sub.o] cos [A.sub.oi] - [[xi].sub.o] sin [A.sub.oi]/tg [z.sub.oi], (2)

[[DELTA].sub.A] = [[xi].sub.o] cos [A.sub.oi] - [[eta].sub.o] sin [A.sub.oi]. (3)

For the reduction of observation data from astronomical to geodetic horizontal coordinate system, the values of azimuth and zenith distance in a formula (1) should be corrected by corrections (2) and (3).

Values of vertical deflection are relatively small. For example, they are below [10.sup.//] in the territory of Lithuania. Therefore, the difference between astronomical and geodetic horizontal coordinate systems can be often ignored when high accuracy is not needed. For example, at zenith angle of 45[degrees], an error in 30 km, distance will be below 1 m due to difference between geodetic and astronomical horizontal coordinates.

3. Space Cartesian geocentric equatorial coordinates of the object

3D Cartesian geocentric equatorial coordinates are tied to the geodetic coordinate system. The origin of the coordinates coincides with the ellipsoidal centre of the geodetic coordinates (Fig), planes XOY, XOZ and YOZ coincide with the planes of the ellipsoidal equator, geodetic prime meridian (L = [0.sup.[omicron]]), and geodetic longitude (L = [90.sup.[omicron]]).

Instrument location at point [P.sub.o] geodetic coordinates are [B.sub.o], [L.sub.o], and [H.sub.o.sup.(g)]. The 3D rectangular coordinates of this point will be (Jakucionis et al. 1999; Zakarevicius 1998):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (4)

where curvature radius of prime vertical

[N.sub.o] = a/[(1 - [e.sup.2] [sin.sup.2] [B.sub.o]).sup.1/2], (5)

a - ellipsoid semi major axis and [e.sup.2] - ellipsoid eccentricity squared.

The horizontal coordinates of the object to be determined are transformed into the 3D Cartesian geocentric equatorial coordinate system (Jakucionis et al. 1999; Zakarevicius 1998):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

The matrix of rotation is (Jakucionis et al. 1999; Zakarevicius 1998):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

The object's 3D location determined in the horizontal coordinate system is tied to the applied Cartesian geocentric equatorial coordinates of the geodetic coordinate system.

4. Determination of geodetic coordinates of object

If the Cartesian geocentric coordinates [X.sub.i], [Y.sub.i], [Z.sub.i] of an object are known, it is possible to compute its geodetic coordinates [B.sub.i], [L.sub.i], [H.sub.i.sup.(g)] (Moritz 1990; Zakarevicius 1998).

Geodetic longitude

[L.sub.i] = arctan [Y.sub.i]/[X.sub.i]. (8)

Geodetic latitude is computed by a method of approximations. First approximation is

[B.sup.(1).sub.i] = arctan [Z.sub.i] cos [L.sub.i]/[X.sub.i] (1 - [e.sub.2]) = arctan [Z.sub.i] sin [L.sub.i]/[Y.sub.i] (1 - [e.sup.2]). (9)

Other approximations:

[B.sup.(k).sub.i] = arctan [T.sub.k-1]/D, (10)

where

D = [X.sub.i]/cos [L.sub.i] = [Y.sub.i]/sin [L.sub.i], (11)

[T.sub.k-1] = [Z.sub.i] - [N.sup.(k-1).sub.i] [e.sup.2] s in [B.sup.(k-1).sub.i], (12)

where [N.sup.(k-1).sub.i]--curvature radius of prime vertical, computed for geodetic latitude [B.sup.(k-1).sub.i].

Approximations are computed until two neighbouring results differ less than the allowable error of determining geodetic coordinates.

Object height in normal direction above the ellipsoid describing geodetic coordinate system:

[H.sup.(g).sub.i)] = [X.sub.i]/cos [L.sub.i] cos [B.sub.i] - [N.sub.i] = [Y.sub.i]/sin [L.sub.i] cos [B.sub.i] - [N.sub.i], (13)

where [N.sub.i]--curvature radius of prime vertical (5), computed for geodetic latitude [B.sub.i].

It is often important to know the normal height of the object (height from main level surface) but not from ellipsoid (13). Normal height:

[H.sup.(n).sub.i] = [H.sup.(g).sub.i] - [h.sub.i], (14)

where [h.sub.i]--geoid height above ellipsoid, which is a function of geodetic coordinates

[h.sub.i] = [F.sub.h]([B.sub.i], [L.sub.i]). (15)

The geoid digital model (15) for the territory of Lithuania is determined with errors below 5 cm (Petroskevicius 2004).

It is often important to know object height above the earth's topographic surface. Object height above the earth's topographic surface will be

[H.sub.i]_= [H.sup.(n).sub.i] - [H.sup.(t).sub.i], (16)

where [H.sup.(t).sub.i]--height of the Earth's topographic surface.

When the digital terrain map is available, the height of the topographic surface is computed as a function of geodetic coordinates

[H.sup.(t).sub.i] = [F.sub.t] ([B.sub.i], [L.sub.i]). (17)

Accuracy determining the height of an object above the Earth's surface is mostly influenced by the accuracy of map digital terrain model.

5. Determination of coordinates of observation instrument

Let us to extend the application of the technological scheme for the determination of geographical coordinates for the solution of the inverse problem, i.e. when 3D polar measurements are produced to the point with the well-known geodetic coordinates and the coordinates of the observation instrument need to be determined (for example, to determine the position of an observation instrument in an aircraft by measurements to a point on the earth with well-known coordinates). Let us assume that the geodetic coordinates of point [P.sub.i] (Fig) are well-known and the coordinates of the observation instrument at point [P.sub.o] need to be determined.

In this case, after measurements from a point with well-known coordinates, one can calculate the conditional intermediate horizontal coordinates (1) of the point on the earth to which measurements are produced. These coordinates will be in the horizontal coordinate system related to point [P.sub.o] because the 3D position of the origin of the coordinate system is unknown.

Furthermore, when one knows the conditional intermediate horizontal coordinates of this point, polar measurement data are deduced into the horizontal coordinate system of point [P.sub.i] with the well-known geodetic coordinates. It is necessary to know distance in space [s.sub.io], zenith angle [z.sub.io] and direction [A.sub.io] of the origin (point [P.sub.i]) of the horizontal coordinate system.

Range is measured directly. It will be equal to the distance measured from point [P.sub.o], i.e., [s.sub.io] = [s.sub.oi]. Direct and inverse directions of zenith distances will not coincide because of the non-coincidence of directions of vertical at points [P.sub.o] and [P.sub.i.] For this reason, this non-coincidence must be taken into account to deduce zenith distance.

Let us assume that ranges are not long (e.g. 20-30 km) and for the simplification of the problem the surface of the ellipsoid may be replaced by the surface of a sphere with the curvature radius R = [square root of M x N], where M and N are the radius of curvature of meridian and radius of the curvature of prime vertical respectively calculated with respect to the geodetic coordinates of point [P.sub.i].

Zenith distance will be

[z.sub.io] = [180.sup.[omicron]] - [z.sub.oi] + [DELTA]z. (18)

Correction [DELTA]z may be calculated by using conditional intermediate coordinates of the horizontal system with origin at point [P.sub.o].

[DELTA]z = arctg ([square root of [x.sujp.2.sub.i] + [y.sup.2.sub.i]]/R + [H.sub.o] + [z.sub.i]), (19)

where R--average radius of the ellipsoid's curvature, H--height of point [P.sub.o]. This height may be measured, in the air vehicle by radar altimeter or even by barometric altimeter. It is acceptable if this height does not have a precise value because the error of the height has little influence on the correction [DELTA]z, e.g. an error having a magnitude of 1 km will give an error of approximately 0,1" in the value of the correction for zenith angle.

In calculations of the value of inverse direction, it is necessary to take into account that the directions of abscissa axes of horizontal coordinates with origins in points [P.sub.o] and [P.sub.i] do not coincide. In accord with the earlier mentioned assumption that distances between these points are not more than several dozen kilometres, the angle between abscissa axes can be found from the formula

[DELTA]A = [arctg ([y.sub.i]/R + [H.sub.o] + [z.sub.i])]tg [B.sub.i]. (20)

Accordingly, the measured direction of point [P.sub.i] in the horizontal coordinate system will be

[A.sub.oi] = [A.sub.oi] [+ or -] [180.sup.]omicron]] + [DELTA]A. (21)

When one has reduced measurement data into the horizontal coordinate system with origin at point [P.sub.i] with well-known coordinates, formulas (1), (4), (5), (6), and (7) can be applied, and after counterchanging indexes 0 and i near the symbols in these formulas, the 3D Cartesian coordinates of the location of the observation instrument can be determined. The application of the scheme proposed here for the reduction of measurement data from one horizontal system to another system with unknown location of its origin will be tied with minor methodical errors. These errors will appear due to the assumption that for moderate distances between points the surface of an ellipsoid can be replaced by the surface of a sphere. For example, for distances up to 30 km, the methodical errors for the 3D Cartesian coordinates of the observation instrument will be less than 2 m.

6. Conclusions

1. Applying intermediate 3D Cartesian horizontal and geocentric equatorial coordinate systems to determine an object's 3D position in the Earth's unified geodetic coordinate system by polar observations is suggested.

2. Using this intermediate rectilinear 3D Cartesian coordinate system, it is possible to simplify the geometric model of polar measurements, avoiding the determination of an object's geometric position linked to the Earth's ellipsoid variable curvature surface.

3. An algorithm to transform polar observations from an astronomic to a geodetic horizontal coordinate system is presented. Angular observation data corrections were computed by making reductions of observation data.

4. There is a presentation of the technological sequence of calculations and algorithms for calculation of 3D Cartesian coordinates of an observation instrument with an initially unknown location (e.g., in an aircraft) when measurements are produced to the point with well-known coordinates.

5. By applying the suggested polar observation method and technologic sequence, as well as algorithms of observation data processing, it is possible to combine observation instrument data with GPS observation data.

DOI: 10.3846/aviation.2010.07

References

Boucher, C.; Altamimi, Z. 2001. ITRS, PZ-90 and WGS 84: current realizations and the related transformation parameters, Journal of Geodesy 75(11): 613-619.

Jakucionis, A.; Zakarevicius, A. 1999. GPS panaudojimo ypatumai tupdymo pagal prietaisus sistemu charakteristiku matavimuose [Peculiarities of GPS usage in ils flight testing], Aviacija [Aviation] 4: 32-34.

Logan, S. A.; Leahy, F. J.; Kealy, A. 2003. Integration of GPS carrier phase and other measurements for kinematic mapping, Journal of Geodesy 76(9-10): 543-556.

Moritz, H. 1990. The figure of the Earth. Theoretical geodesy and the Earth's Interior. Karlsruhe: Wichmann. 279 p.

Petroskevicius, P. 2004. Gravitacijos lauko poveikis geodeziniams matavimams [The effect of the gravitational field on geodetic observations]. Vilnius: Technika. 292 p. ISBN 9986-05-738-8.

Skeivalas, J. 2008. GPS tinklu teorija ir praktika [Theory and practice of GPS networks]. Vilnius: Technika. 288 p. ISBN 978-9955-28-228-0.

Zakarevicius, A. 1994. Dabartiniu vertikaliu Zemes plutos judesiu Lietuvos teritorijoje tyrimas [Investigation of present vertical crust movements of the Earth's in the territory of Lithuania]. Vilnius: Technika. 276 p. ISBN 9986-05-097-9.

Zakarevicius, A. 1998. Navigaciniu duomenu perskaiciavimas i WGS-84 koordinaciu sistema [Navigation data conversion to coordinate system WGS-84], Aviacija [Aviation] 3: 100-102.

Zakarevicius, A. 2003. Dabartiniu geodinaminiu procesu Lietuvos teritorijoje tyrimas [Investigation of recent geodynamic processes in the territory of Lithuania]. Vilnius: Technika. 195 p. ISBN 9986-05-691-8.

Algimantas Zakarevicius (1), Vladislovas Ceslovas Aksamitauskas (1), Algimantas Jakucionis (2), Arminas Stanionis (1)

(1) Vilnius Gediminas Technical University, Sauletekio al. 11, LT-10223 Vilnius, Lithuania (2) Antanas Gustaitis Aviation Institute of Vilnius Gediminas Technical University, Rodunes kelias 30, LT-02187 Vilnius, Lithuania E-mail: (1) gkk@vgtu.lt; (2) Algimantas.Jakucionis@vgtu.lt

Received 3 May 2010; accepted 20 May 2010

Algimantas ZAKAREVICIUS, Prof Dr Habil

Date and place of birth: 27 April 1939 in Lithuania.

Education: Kaunas Polytechnic Institute (now Kaunas University of Technology); geodetic engineer (1965).

Affiliations and functions: 1973--doctoral degree at Vilnius University, 2000--Dr Habil degree at VGTU.

Research interests: investigation of the recent geodynamic processes, formation of geodetic networks.

Publications: 170, including 3 monographs.

Present position: professor, Department of Geodesy and Cadastre, Vilnius Gediminas Technical University, Sauletekio al. 11, LT-10223 Vilnius, Lithuania, +37052370630, Algimantas.Zakarevicius@vgtu.lt.

Ceslovas AKSAMITAUSKAS, Assoc Prof Dr

Date and place of birth: 2 July 1947 in Lithuania.

Education: Vilnius Civil Engineering Institute; geodetic engineer (1975).

Affiliations and functions: 1982--doctoral degree at Lviv Polytechnic Institute.

Research interests: geodetic instruments, investigations of deformations, GPS.

Publications: 80.

Present position: associate professor, head of Department of Geodesy and Cadastre, Vilnius Gediminas Technical University, Sauletekio al. 11, LT-10223 Vilnius, Lithuania, +37052744701, caks@vgtu.lt.

Algimantas JAKUCIONIS Assoc Prof Dr

Date and place of birth: 6 February 1942 in Lithuania.

Education: Kaunas Polytechnic Institute (now Kaunas Technological University); radio electronics engineer (1964).

Affiliations and functions: 1974--doctoral degree in technical sciences at Kaunas Politechnic Institute.

Research interests: GNSS applications, measurement and data interpretation methods for automated flight inspection of NAVAIDS.

Publications: 40.

Present position: associate professor, head of Department of Avionics, Antanas Gustaitis Aviation Institute of Vilnius Gediminas Technical University, Rodunes kelias 30, LT-02187 Vilnius, Lithuania, +37052744817, Algimantas.Jakucionis@vgtu.lt.

Arminas STANIONIS, Assoc Prof Dr

Date and place of birth: 15 April 1976 in Lithuania.

Education: Vilnius Gediminas Technical University (VGTU); bachelor's degree in Geodesy (2000), master of science in Measurement Engineering (2002).

Affiliations and functions: 2005--doctoral degree in Measurement Engineering. Research interests: investigation of geodynamic processes, GIS, investigations of deformations.

Publications: 30.

Present position: associate professor, Department of Geodesy and Cadastre, Vilnius Gediminas Technical University, Sauletekio al. 11, LT-10223 Vilnius, Lithuania, +37052744703, Arminas.Stanionis@vgtu.lt.
COPYRIGHT 2010 Vilnius Gediminas Technical University
No portion of this article can be reproduced without the express written permission from the copyright holder.