CC BY
10
УДК 528.21
РАСЧЕТ МОДЕЛИ ЛОКАЛЬНОГО КВАЗИГЕОИДА (УЛАН-БАТОР),
ОЦЕНКА ПАРАМЕТРОВ И ПОДХОДЫ В ОПРЕДЕЛЕНИИ ГРАВИТАЦИОННОГО
ПОЛЯ
Катерина Морозова
Университет Латвии, Институт геодезии и геоинформатики, LV-1586, Латвия, г. Рига, Бульвар Раиниса, 19; Рижский технический университет, LV-1048, Латвия, г. Рига, ул. Межа, 1, магистр инженерных наук, аспирант, исследователь, тел. (371)26-812-422, e-mail: katerina.morozova@rtu.lv
Райнер Ягер
Высшая школа Карслруе - Университет прикладных наук, Институт прикладных исследований, 76133, Германия, г. Карлсруе, ул. Мольткештрассе, 30, доктор инженерных наук, профессор, глава лаборатории ГНСС и навигации, тел. 49(0)152-533-103-28, e-mail: reiner.jaeger@web.de
Саандаар Миджидордж
MonMap Инженерный сервис Co., 210646, Монголия, г. Уланбатор, ул. Сеул, 12/6, доктор инженерных наук, генеральный директор, тел. (976)99-11-32-61, e-mail: msaandar@monmap.mn
Гунарс Силабриедис
Университет Латвии, Институт геодезии и геоинформатики, LV-1586, Латвия, г. Рига, Бульвар Раиниса, 19, доктор инженерных наук, директор института, тел. (371)27-476-220, e-mail: gunars.silabriedis@lu.lv
Янис Балодис
Университет Латвии, Институт геодезии и геоинформатики, LV-1586, Латвия, г. Рига, Бульвар Раиниса, 19, доктор физических наук, ведущий исследователь, тел. (371)29-298-150, e-mail: Janis.balodis@lu.lv
Янис Каминскис
Рижский технический университет, LV-1048, Латвия, г. Рига, ул. Межа, 1, доктор инженерных наук, руководитель кафедры; Университет Латвии, Институт геодезии и геоинформатики, LV-1586, Латвия, г. Рига, Бульвар Раиниса, 19, ведущий исследователь, тел. (371)27-476-220, e-mail: Janis.kaminskis@lu.lv
В работе описан метод и результаты расчета локальной модели квазигеоида для региона Уланбатора, UBQGEOID2018, основанный на глобальной гравитационной модели EGM2008 [1], ГНСС и опорных точках нивелирования в Балтийской системе высот 1977, а также производных данных линий отвеса. Для расчетов использовалось програмное обеспечение DFHRS в. 4.4, которое позволяет напрямую конвертировать эллипсоидальные высоты в нормальные высоты, основываясь на параметрическом моделировании опорной поверхности высот. В статье приводится принцип измерений отвесных линий зенитной камерой на основе матрица ротаций: матрица ротации между локальной астрономической вертикальной системой и локальной геодезической системой также представлены. Описаны следующий этап развития ПО и концептуальные формулы.
Ключевые слова: линия отвеса, DFHRS, EGM2008, геодезия, геоид, квазигеоид, геофизика, модели геопотенциала, гравитационное поле, зенитная камера.
ULAANBAATAR QGEOID COMPUTATION, PARAMETER ESTIMATION AND OPTIMIZATION CONCEPTS FOR GRAVITY FIELD DETERMINATION
Katerina Morozova
University of Latvia, Institute of Geodesy and Geoinformatics, 19, Boulevard of Rainis St., of. 408, Riga, LV-1586, Latvia; Riga Technical University, 1, Meza St., of. 116, Riga, LV-1048, Latvia, MSc, Ph. D. Student, Researcher, phone: (371)26-812-422, e-mail: Katerina.morozova@rtu.lv
Reiner Jäger
Hochschule Karlsruhe - University of Applied Sciences, Institute of Applied Research, 30, Moltkestrasse, of. 212, Karlsruhe, 76133, Germany, D. Sc., Professor, Head of Laboratory GNSS and Navigation, phone: (490)152-533-103-28, e-mail: reiner.jaeger@web.de
Saandaar Mijiddorj
MonMap Engineering Services Co., Ltd, 12/6, Seoul St., of. 502, Ulaanbaatar, 210646, Mongolia, D. Sc., Director General, phone: (976)99-11-32-61, e-mail: msaandar@monmap.mn
Gunars Silabriedis
University of Latvia, Institute of Geodesy and Geoinformatics, 19, Boulevard of Rainis, of. 406, Riga, LV-1586, Latvia, Dr. Sc. ing., director, phone: +371 27 476 220, e-mail: gunars.silabriedis@lu.lv
Janis Balodis
University of Latvia, Institute of Geodesy and Geoinformatics, 19, Boulevard of Rainis St., of. 406, Riga, LV-1586, Latvia, D. Sc., Leading Researcher, phone: (371)29-298-150, e-mail: Janis.balodis@lu.lv
Janis Kaminskis
Riga Technical University, 1, Meza St., of. 116, Riga, LV-1048, Latvia, D. Sc., Head of Department; University of Latvia, Institute of Geodesy and Geoinformatics, 19, Boulevard of Rainis St., of. 407, Riga, LV-1586, Latvia, Leading Researcher, phone: (371)27-476-220, e-mail: Janis.kaminskis@lu.lv
The article describes the method and computation results of Ulaanbaatar region quasi-geoid model, UBQGEOID2018, based on global gravity field model EGM2008 [1], GNSS and levelling points in Baltic Height system 1977 and derived deflections of vertical data. The DFHRS (Digital Finite-element Height Reference Surface) software v.4.4 [2] has been used for this purpose, which allow the direct conversion of ellipsoidal heights to normal heights and based on parametric modelling of the HRS. The principle of vertical deflections measurements by digital zenith camera is included in this paper based on rotation matrices: the rotation matrix between local astronomical vertical system and local geodetic vertical system is introduced. The next stage of the development version 5.x and related research is also described and the conceptual formulas are introduced.
Key words: deflection of vertical, DFHRS, EGM2008, geodesy, Geoid, geophysics, geopotential models, gravity field, Quasi-Geoid, zenith camera.
Introduction
In the era of modern technologies and GNSS developments the precise quasi-geoid model is necessary for different engineering needs, as it allows the determina-
tion of normal height much faster in comparison to levelling and directly from GNSS. This article describes the software for quasi-geoid determination based on parametric modelling, as well as further version based on Spherical-Cap-Harmonics (SCH) modelling. The example of quasi-geoid model for Ulaanbaatar region and computation results are introduced. The theory of deflections of vertical measurements by digital zenith camera is also included.
1. Principle of the DFHRS software
The principle of a GNSS-based height determination H, requires submitting the GNSS-height h to the DFHRS (B, L, h)-correction, in order to receive physical height H and it reads:
H = h - DFHRS(p, Am | B, L, h) = h - (NFEM(p | B, L) + Am ■ h (2-1)
The DFHRS-correction DFHRS (B, L, h) is provided by means of a DFHRS database (DFHRS_DB), which contains the HRS polynomial parameters and the scale difference (p, Am) together with the mesh-design information. The mathematical model for observation groups in a common least squares computation (GauB-Markov-Model) for the evaluation of the DFHRS_DB parameters p and Am is given by formulas (2-2a-f) [2, 3].
Functional Model
h + v = H + h ■Am + NFEM (p | x, y), with NFEM(p | x, y) =: f (x, y) ■ p
Observation Types and Stochastic Models Uncorrelated ellipsoidal height h observations. Covariance
matrix C
h
2
diag
i
(2-2a)
Ng (B, L)J + v = f(x, y)T ■ p + ÔNg (dJ )
T
J
Correlated geoid height observations. With a given real co-
variance matrix CN or CN
nG nG
evaluated from a synthetic co-variance function.
(2-2b)
tJ + v = -fBT/(M (B) + h) ■ p +
y + v = - fT / (N(B) + h) ■ cos(B)) ■ p + +dn(d-Jz )
Observations of deflections from the vertical £). Pair-wise correlated or uncorrelated (2-2c) in case of astronomical observations. Correlated if derived from a gravity potential model.
(2-2d)
Uncorrelated standard height H
TT TT observations with covariance _ N
H+v = H 1 (2-2e)
matrix C^ = diag(ct^ )
i
Continuity condition equations (1d) introduced as uncorrelated C + v = C(p) so-called pseudo observations (2-2f)
with accordingly small variances and high weights.
With ffi and fj^ we introduce the partial derivatives off(x(B, L), y(B, L)) (2-2c)
with respect to the geographical coordinates B and L. M(B) and N(B) mean the radius of meridian and normal curvature at a latitude B. The continuity of the resulting HRS
representation NFEM(p | x, y) = f (x, y) x p over the meshes (fig. 1, thin blue
lines) is automatically provided by the continuity equations C ( p ) (2-2f). A number
of identical fitting-points (B, L, h; H) are introduced by the observation equations (2-2a) and (2-2e) (fig. 1, green triangles). In the practice of DFHRS_DB evaluation, one or a number of different geoid-/GPM such as the EGG97 or EGM 2008 are used in a least squares estimation related to the mathematical model (2-2a-f), which is implemented in the DFHRS-software 4.4. To reduce the effect of medium- or longwave systematic shape deflections, namely the natural and stochastic "weak-shapes", in the observations N and (£, from geoid- or GPM models, these observations are subdivided into a number of patches (fig. 1, thick blue lines). These patches are related to a set of individual parameters, which are introduced by the datum
parametrizations BNq ( dj ) (2-2b) and d£(d * );^ ) (2-2c, d). In this way,
it is of course possible to introduce geoid height observations and vertical deflections from any number of different geoid- or GPM models in the same area, or observed vertical deflections [2, 3].
2. Computation results of DFHRS-based Ulaanbaatar Region Quasi-Geoid for the Baltics Height System
In order to compute the DFHRS_DB for Ulaanbaatar 94 Identical points (ellipsoidal heights h and normal heights H in Baltic Height system) together with the EGM2008 geopotential model data were used. EGM2008 is a spherical harmonic model of the earth's external gravitational potential in degree and order of 2160, with additional spherical harmonic coefficients extending up to degree of 2190 and order of 2160 that offers a spatial resolution of 9 km. EGM2008 incorporates improved 5 x 5 min gravity anomalies, altimetry-derived gravity anomalies and has benefited from the latest GRACE based satellite solutions [4].
Fig. 1. Computation design of DFHRS (meshes - thin blue lines, patches - thick blue lines, fitting points - green triangles)
For meshing the area, mesh size of 5 x 5 km was chosen (fig. 1, thin blue lines). Total amount of meshes - 1536. The total number of patches is 5 (fig. 1). One patch must contain at least 4 fitting points. As points of the region are not homogenously located, patches, were not introduced in approximately the same size, but according to the location of the points. As geoid datum 3 translations and 3 rotations were introduced, additionally derived deflections of the vertical from the EGM2008 model were used (see fig. 2).
The identical points and the EGM2008 geoid undulations were introduced together with the continuity conditions into a least squares estimation of the so-called "DFHRS production". The calculation has been done using the DFHRS v. 4.4. software. 88 normal height points H of the Baltic heights system could be used and were confirmed in the statistical testing (data-snooping) with the assumed standard deviation of 1 cm. 6 points - 4039, 216, 230, 5051, 509 and 22 were excluded from the computations because of gross errors. For 4 points (270, 1710, 1757 and GR70/70) the normal heights H were changed in comparison to the previous data package provided in 2017, the normal height for one point (1682) was used from previous data package.
Fig. 2. DFHRS-software 4.4 computation dialog
The partial adjustment protocol of the DFHRS-software 4.4. with the observation residuals, statistical testing of the height fitting points is depicted in table 1.
Table 1
Final DFHRS software adjustment protocol
Characteristics:
EV: Redudancy factor
NV: Normalized residuals, test size a priori
t post: Test size a posteriori
GF: Estimated gross error is issued in case of exceeding
the critical value by nv, bzw. T_post.
Probability of error Alpha: 5 %
Critical value a priori: 3.841552 degrees of freedom: infinity
Critical value a posteriori: 3.841549 degrees of freedom: 102887
Point number Height/Target sys. Res. EV NV t_post REPRO
[m] [m] [%]
137 1168. 601 0. 00009 18. 88 0. 0 0. 5 -0. 001
253 1441. 619 0. 00449 22. 69 1. 3 22. 2 -0. , 020
268 1459. 774 -0. 00176 21. 66 0. 5 8. 9 0. 008
282 1318. 717 -0. 00275 19. 97 0. 9 14. 5 0. 014
505 1433. 682 -0. 00366 22. 85 1. 1 18. 0 0. 016
The rest of the table 1
Point number Height/Target sys. Res. EV NV t_post REPRO
[m] [m] [%]
1598 1418. 122 0. 00106 16. 46 0. 4 6. 1 -0. 006
1710 1478. 901 0. 00345 22. 20 1. 0 17. 3 -0. 016
1731 1608. 961 0. 00347 20. 96 1. 1 17. 9 -0. 017
1747 1313. 276 -0. 00114 23. 60 0. 3 5. 5 0. 005
1757 1220. 986 0. 00231 20. 85 0. 7 11. 9 -0. 011
2324 1263. 574 -0. 00330 23. 20 1. 0 16. 1 0. 014
2329 1246. 173 -0. 00198 23. 35 0. 6 9. 6 0. 008
4750 1456. 719 -0. 00098 21. 87 0. 3 4. 9 0. 004
5006 1419. 172 0. 00655 23. 67 1. 9 31. 8 -0. , 028
5019 1425. 403 0. 00617 23. 32 1. 8 30. 2 -0. , 026
Excluded fitting points from the computations are depicted in table 2.
Table 2
Eliminated Error Points
Probability of error Alpha: 5 %
Critical value a priori: 3.841552 degrees of freedom: infinity
Critical value a posteriori: 3.841549 degrees of freedom: 102889
Point number Height/Target sys. Res. EV NV t_post REPRO
[m] [m] [%]
509 1443.301 -0.05132 27.61 14.0* 180.2 0.186**
!!! ---> gf: 0.186 m <--- !!!
5051 1354.620 -0.03123 27.97 8.4* 99.4 0.112**
!!! ---> gf: 0.112 m <--- !!!
22 1228.700 -0.01362 13.73 5.3* 79.0 0.099**
!!! ---> gf: 0.099 m <--- !!!
4039 1435.117 -0.01634 22.59 4.9* 73.6 0.072**
!!! ---> gf: 0.072 m <--- !!!
216 1353.229 0.01688 26.19 4.7* 70.5 -0.064**
!!! ---> gf: -0.064 m <--- !!!
230 1264.222 0.01978 25.68 5.6* 84.2 -0.077**
!!! ---> gf: -0.077 m <--- !!!
3. Conclusions and results for computed Ulaanbaatar QGeoid model
The present DFHRS was calculated on the basis of the EGM2008 geoid and 88 identical reference points. The accuracy of the identical points was confirmed with 1.0 cm, so the QGeoid of the Ulaanbaatar region has an estimated 1 -3 cm accuracy within the area of the outer ring polygon-line of the fitting-points. The DFHRS_DB can be used by the software DFHBF-Tools to compute the QGeoid-height N, and so
the Normal Heights H from the input of a 3D GNSS-position (B, L, h) or (X, Y, Z), and in order to set up a respective QGeoid 2018 grid for the Baltic Height System in the Ulaanbaatar Region. Especially for the borders of the Region (fig. 1) additional vertical deflection observations made by digital zenith camera [5, 6] are recommended. In that way, the 1-3 cm accuracy will hold for the whole area (fig. 1).
4. Zenith camera and determination of deflections of the vertical
The basic component are imaging sensors (CCD cameras) to track celestial objects or stars, respectively. If we suppose, that the imaging sensor system and the respective platform p is already aligned or identical with the body system b, we have for p and the direction vector S of the body system [7]:
p =
SI
x xq
y - yQ f
and 4 =p = ^ •
p p
x xq
y - yQ f
(5-1)
rLav (O, A) =
with (x, y)t UTc - Observed image coordinates of a star S(5, a) at time t_UTC;
(xQ, yQ) - Principal point of the image, f- Focal length of the sensor (see, fig. 3). Further we have
- cos A- sin O - sin A- sin O cos O
- sin A + cos A Q . (5-2)
cos A- cos O sin A- cos O sin O
The astronomical position is described with (O, A) and the geographical GNSS-position with (B, L) leading to
- cos L • sin B - sinL • sin B cos B
- sin L cos L Q . (5-3)
cos L • cos B sinL • cos B sinB
R LGV (B, L) =
From (5-2) and (5-3) we get
-nLAV olav lgvt
R LGV = R e • (Re ) :
(5-4)
with
R LGVV=R LGVV (B, L,^)=
^sinB sinO cos( A - L) + cos B cos O - sin O sin(A- L)
sin B sin( A - L) cos B sin O - sin $cos O cos( A - L) cos(A - L) + cos O sin(A- L)
sin B cos O - cos B sin O cos( A - L) - cos B sin( A - L) cos B cos O cos( A - L) + sin B sin O
\T
(5-5)
Fig. 3. Modern Star tracker CT-602 as produced by Ball Aerospace's CT-602
e s
With the star coordinates r ' at the observation time t_UTC we have
rLGV = rLgv (B, L) -re's.
(5-6)
All in all the general model for the vertical surface deflections determination the equation reads:
rLGV (5 - 6) - rLGV (5 - 8«, b) = 0
(5-7)
with
,.LGV nLAV
RLGV (b, l^)t • rLav (r = o, p = o, y) • r°SI = o.
SI
(5-8a)
The matrix
R
LAV
^cos p cosy sin r sin p cos y - cos r sin y cos r sin p cos y + sin r siny^
cos p sin y sin r sin p siny + cos r cos y cos r sin p siny - sin r cos y v - sin p sin r cos p cos r cos p J
(5-8b)
is by the horizontation (r = 0, p = 0) of the zenith camera platform in the local LAV using an inclinometer sensor to R^ (r = 0, p = 0, y). The heading y is approximately
LGV
known, but remains an unknown of the parameter estimation. For RMv in (5-8a) we can also use [7]
b
R
LGV LAV
R
LGV LAV
1 tan B
tan B 1
(5-9)
5. Next stage of the software development - DFHRS v. 5.x
The extension of DFHRS concept and software to physical observation types -such as terrestrial, air- or space-borne gravity measurements or physical observation types taken from geopotential models, e. g. EGM 2008 - is based on a regional adjusted spherical cap harmonic parameterization (ASCH) of the Earth's gravitational potential (V) [2, 8, 9]:
V ( r,X ', 0 ')= I
k=0
k max/a\n(k)+1 k
V r y
I (Cn(k)m C0S'+ Sn(k)m sin ')Pn(k),m (c0S 0') (6-1)
m=0
New adjustment-based approach enables estimation of coefficients (C
( k )
n( k )m
S.
iir\m for regional ASCH model V as functions of coefficients (C , S ) of a
..( k )m » v n,m n,m
global geopotential model. The estimated coefficients (Cn^k)m, Sn^k)m ) can be introduced as so-called direct observations in the integrated approach, and thus we have:
C '
n(k ),m
(t ) + v = C '
n (k), m
and S 'n(k ),m (t ) + V = S
n(k),m .
(6-2)
In the so-called integrated DFHRS approach we have the following observation equation for a gravity observation:
lgv = GM_ «
^av r2 k=oV r
^ay(k )+1 k
(n(k) + 1) I (C'n(k)),m* C0S'+ S 'n(k),m" sin0 ' Pn(k),m(c0s0 0 . (6-3)
m=0
By introducing the disturbance potential applied to the Bruns theorem and Molodenski's theory, we obtain the observation equation for fitting-points (h - H) converted to quasi-geoid heights NOG and vertical deflections (£, n)p at measured at the earth surface by zenith camera (fig. 1) at a point P reading [2], [9], [10]:
QN0G QB S p =--0G. — + dNCurv
QB Qs
N
h - H = Nqg = ^ Y0
-—( TP ) — + dNCurv
QB Yq QsN
-1
Yq ■ (M + h) QB
(6-4)
QT
(—) P + dNCurv.
(6-5)
p qnqg QL Q 1 QL "H P = —■ — = tt (— Tp )
( QT ) p
QL QsE QL Yq QsE Yq ■ (N + h) ■ cos B QL
(6-6)
One further research topic in the DFHRS-project will be dealing with the optimal design (1st Order Design) of the observation type of gravity observations (6-3) and vertical deflection observations (6-5,6).
Conclusions
The quasi-geoid model for Ulaanbaatar region has been computed. The accuracy of the model is evaluated by 1-3 cm. As levelling data are not homogeneously provided in the region of interest, it would be necessary to use digital zenith camera for vertical deflection determination for quasi-geoid improvement, as well as it allows additional check of normal heights. ASCH modelling in terms of integrated geodesy allow the combination of both geometrical and physical data, moreover this method is much faster in comparison to SH. Implementation of vertical deflections observations in terms of ASCH gives additional improvement of quasi-geoid and gravity field determination.
REFERENCES
1. International Center for Global Gravity Field Models (2018) http://icgem.gfz-potsdam.de/ICGEM/
2. DFHBF-Website (2000-2018): www.dfhbf.de
3. Jäger R., Kaminskis J., Strauhmanis J. and Younis G. (2012): Determination of quasi-geoid as height component of the geodetic infrastructure for GNSS positioning services in the Baltic States," Latvian J. of Physics and Technical Sciences 3, pp. 5-15.
4. Pavlis N. K, Holmes S. A., Kenyon S. C, Factor J. K. (2008). An Earth Gravitational model to degree 2160: EGM2008, General Assembly of the European Geosciences Union, Vienna, Austria
5. A. Zarins, A. Rubans, and G. Silabriedis, (2016) Digital zenith camera of the University of Latvia, Geodesy and Cartography, 42:4, pp. 129-135. http://dx.doi.org/10.3846/20296991. 2016.1268434.
6. Morozova, K., Balodis, J., Jäger, R., Zarins, A., Rubans, A. Digital Zenith Camera's Results and Its Use in DFHRS v.4.3 Software for Quasi-geoid Determination (2017). From: 2017 Baltic Geodetic Congress (BGC Geomatics), Polija, Gdansk, 22.-25. June, 2017. Piscataway: IEEE, 2017, 174.-178.lpp. ISBN 978-1-5090-6041-2. e-ISBN 978-1-5090-6040-5. Available from: doi:10.1109/BGC.Geomatics.2017.74
7. Jekeli, C (2000): Inertial Navigation Systems with Geodetic Applications. De Gruyter.
8. Younis, G. (2013): Regional Gravity Field Modeling with Adjusted Spherical Cap Harmonics in an Integrated Approach. Schriftenreihe Fachrichtung Geodäsie der Technischen Universität Darmstadt (39). Darmstadt. ISBN978-3-935631-28-0 (2013)
9. G. K. A. Younis, R. Jäger, and M. Becker, (2011) Transformation of global spherical harmonic models of the gravity field to a local adjusted spherical cap harmonic model, Arabian Journal of Geosciences. DOI 10.1007/s12517-011-0352-1.
10. Morozova, K., Jäger, R., Balodis, J., Kaminskis, J. (2017) Software Development and Its Description for Geoid Determination Based on Spherical-Cap-Harmonics Modelling Using Digital-Zenith Camera and Gravimetric Measurements Hybrid Data. IOP Conference Series: Materials Science and Engineering, 2017, Vol.251, pp.1-10. ISSN 1757-8981. e-ISSN 1757-899X. Available from: doi:10.1088/1757-899X/251/1/01206.
© К. Морозова, Р. Ягер, С. Миджидордж, Г. Силабриедис, Я. Балодис, Я. Каминскис, 2018
