G i I 60/ ^^ V GEODETSKI VESTNIK | letn./ Vol. 60 | št. / No. 2 | 1 določitev primernega definition of i geodetskega datuma appropriate geodetic t z uporabo robustnih datum using robust statističnih metod statistical methods Aleš Marjetic, Klemen Kregar 9 UDK: 528.1:528.3:519.246 Klasifikacija prispevka po COBISS.SI: 1.01 Prispelo: 18. 3. 2016 Sprejeto: 13. 5. 2016 DOI: 10.15292/geodetski-vestnik.2016.02.212-226 SCIENTIFIC ARTICLE Received: 18. 3. 2016 Accepted: 13. 5. 2016 IZVLEČEK ABSTRACT Pravilna določitev geodetskega datuma je nujen pogoj za pravilno določitev premikov točk. Obstajajo različne metode deformacijske analize, katerih bistvo je pravilna določitev tistih točk, ki niso podvržene premikom in ki lahko predstavljajo primerno koordinatno osnovo za izračun realnih vrednosti premikov ostalih točk. Te metode v večini temeljijo na statističnemu testiranju. Predstavljajo celovito mreže in omogočajo določitev statistično značilnih premikov. Če poznamo lastnosti transformacije med rešitvami premikov, ki so določene na podlagi različnih definicij geodetskih datumov, lahko problem definiranja primernega geodetskega datuma rešujemo nekoliko drugače. V obravnavanem članku smo se osredotočili na določitev ustrezne utežne matrike E v transformaciji S. Pri tem smo uporabili splošno znane metode robustne statistike. Robustnost treh izbranih metod smo preverili na dveh različnih situacijah izbranih premikov v obravnavani geodetski mreži ter rezultate primerjali tudi na izbranem primeru geodetske mreže z rezultati klasičnih metod deformacijske analize. The correct determination of geodetic datum is an obligatory condition for the proper determination of the point displacements. There are various methods of deformation analysis focused on the right identification of stable points, which may define an appropriate coordinate basis for calculating the displacements of other points. These methods are largely based on the statistical testing and represent a comprehensive and complex analysis of change of geometry of geodetic network and allow definition of statistically significant displacements. Knowing the characteristics of the transformation between the solutions ofdisplacements that are based on different definitions of geodetic datums, the problem of defining an appropriate geodetic datum can be solved in a slightly different way. In this article we have focused on the problem of determining the appropriate weighting matrix E in the model of S-transformation. We used the generally known methods of robust statistics. The robustness of the three selected methods were tested on two different situations of preselected displacements in the considered geodetic network and the results on selected case of geodetic network with results of conventional methods of deformation analysis were compared. KLJUČNE BESEDE KEY WORDS geodetski datum, transformacija S, robustna statistika, premik, vplivna funkcija, iterativno prilagajanje uteži Aleš Marjetic, Klemen Kregar | DOLOČITEV PRIMERNEGA GEODETSKEGA | 212 | ROBUST STATISTICAL METHODS | 212-226 | geodetic datum, S transformation, robust statistics, displacement, influence function, IWP - iterative weighted projection j STATISTIČNI h DEFIN J OF APPROPRIATE GEODETIC USIN GEODETSKI VESTNIK |60/2| 1 UVOD Deformacije v geodetskem smislu v splošnem obravnavamo kot premike posameznih točk objekta, ki so položajno določene z meritvami (dolžin, kotov, višinskih razlik, GNSS-opazovanj itd.). Meritve ne podajajo vseh informacij o geodetskem datumu. Z izravnavo geodetske mreže po metodi najmanjših kvadratov (v nadaljevanju: MNK) izračunamo rešitev za vektor koordinatnih neznank točk mreže v posamezni epohi t. Če predpostavimo, da so opazovanja med posameznimi izmerami med seboj neodvisna, lahko vektorsko polje premikov predstavimo z razliko koordinat točk mreže x1 in x2 med dvema epohama t in t2: u = xx - x2. (1) Posamezni rezultati izravnave geodetske mreže (x1 in x2) ter u so seveda datumsko pogojeni, kar pomeni, da vsebujejo tudi vse informacije o datumu geodetske mreže v posamezni terminski izmeri. Tako lahko vektor u v začetni fazi imenujmo samo vektor koordinatnih sprememb. Pogoji, da lahko spremembe koordinat definirajo tudi vektorsko polje premikov, so opisani v Xu et al., 2000, in Marjetič, 2013. Najpomembnejša pa je zahteva, da moramo zagotoviti enak in smiselno definiran geodetski datum v obeh časovno ločenih izmerah, ki ju primerjamo. Enak koordinatni sistem v dveh ločenih izmerah bi teoretično lahko zagotovili, če bi imeli enako geometrijo mreže, enak tip in število opazovanj, enak instrumentarij, enake vremenske pogoje, enake vrednosti danih količin ter način obdelave opazovanj (Sterle, 2007). To pa je tako rekoč nemogoče zagotoviti. Zato iz vseh navedenih razlogov vektorja u ne moremo obravnavati kot vektor premikov točk, ampak samo vektor koordinatnih sprememb. Opisane težave po nekaterih standardnih metodah deformacijske analize, kot so Delft, Hannover, Karlsruhe, Muenchen in Fredericton (Chrzanowski et al., 1983; Van Mierlo, 1978; Welsch in Zhang, 1983; Vrečko in Ambrožič, 2013), rešujejo s primerjavo koordinat točk izravnave proste mreže dveh izmer. S tem zagotovijo enakost geodetskih datumov, in sicer z identično datumsko matriko, ki zagotavlja ustrezne notranje datumske vezi v geodetski mreži. S postopki statističnega testiranja se določijo domnevno mirujoče točke geodetske mreže, s katerimi lahko na novo definiramo geodetski datum s transformacijo S (Marjetič in Stopar, 2007). Tako transformiramo vektor koordinatnih sprememb preostalih točk v datum domnevno mirujočih točk. V članku problem definicije primernega geodetskega datuma obravnavamo nekoliko drugače. Prek lastnosti transformacije S in osnov robustne statistike lahko iterativno lociramo množico točk, ki se med dvema terminskima izmerama niso premaknile, in točke, ki ne ohranjajo položaja med izmerama in so z vidika robustne statistike grobo pogrešena opazovanja. 2 DOLOČITEV GEODETSKEGA DATUMA Z ROBUSTNIMI STATISTIČNIMI METODAMI Pri iskanju optimalnega geodetskega datuma z metodami robustne statistike (Chen et al., 1990) obravnavamo funkcionalni model transformacije S. Imamo vektor premikov točk u in pripadajočo matriko kofaktorjev Q^u kot razliko rezultatov izravnave (proste mreže) dveh terminskih izmer. Pogoj obeh rešitev terminskih izmer je, da sta izračunani v istem geodetskem datumu. Na podlagi lastnosti transformacije S (Baarda, 1981; Marjetič in Stopar, 2007) lahko rešitev u transformiramo iz datuma i v poljubno izbran geodetski datum j. Model transformacije S ima obliko: Aleš Marjetic, Klemen Kregar | DOLOČITEV PRIMERNEGA GEODETSKEGA DATUMA Z UPORABO ROBUSTNIH STATISTIČNIH METOD | DEFINITION OF APPROPRIATE GEODETIC DATUM USING ROBUST STATISTICAL METHODS | 212-226 | | 215 | |60/2| GEODETSKI VESTNIK u j = S . u; = ^ I - H ( H T E j H )-1 H TEj Jj u, (2) Q' = S .Q; ST, J uu J ' kjer indeks i in j označuje geodetski datum, ter: I — enotska matrika dimenzije 2 m x 2 m, kjer je m število točk v mreži, HT = 1 0 1 0 .. .1 0 0 1 0 1 .. . 1 0 0 0 0 _ 0 - yi x1 - Jl x2 .. . Jm x Xm 0 0 0 0 0 x1 J1 x2 Jl .. .x . Xm Jm — matrika dimenzije 4 x 2 m, ki izhaja iz pogojev za definiranje notranjih vezi za določitev geodetskega datuma (Caspary, 1988) oziroma matrika Helmertove transformacije za vse točke geodetske mreže, kjer prvi dve vrstici podajata zahtevo, da se mreža v povprečju ne premakne, tretja, da se ne zasuka, in četrta, da ne spremeni merila po izravnavi. E.—matrika dimenzije 2 m x 2 m, katere izvendiagonalni elementi so enaki 0, na diagonali pa so vrednosti 1 samo na mestih, ki pripadajo posamezni koordinatni komponenti, ki je dana količina za definiranje geodetskega datuma i. Če zapišemo Gauss-Markov funkcionalni model geodetske mreže, imamo (Teunissen, 2003): v + BA = f= d - l. (3) V teoriji izravnave po metodi najmanjših kvadratov iščemo rešitve, ki izpolnjujejo pogoj najmanjše vsote kvadratov popravkov opazovanj. Te rešitve dobimo, ko izpolnimo pogoj ortogonalnosti med vektorjem popravkov opazovanj in prostorom matrike B (Teunissen, 2003), pri čemer upoštevamo različen vpliv opazovanj, ki ga podajamo z matriko uteži P: BTPv = 0, (4) BTP(f - BA) = 0, iz katerega izhajajo normalne enačbe: BTPBA = BTPf. (5) Iz (4) in (5) izhajajo rešitve izravnave po metodi najmanjših kvadratov (tu jih posebej označujemo s strešico): A = (BtPBBTPf = N-1t, (6) i = d - BA , v = f - BA = f - B (B tPB)-1 BTPf = f I - B (b tPB)-1 BTP j f . Primerjava modelov (2) in (6) kaže, da je funkcionalni model transformacije S (2) dejansko Gauss-Markov model, ki povezuje »opazovanja« u, »popravke« u. ter »neznanke« vektorja transformacijskih parametrov t5: u + Htc = u, j ^ i kjer je rešitev za t5: Aleš Marjetič, Klemen Kregar | DOLOČITEV PRIMERNEGA GEODETSKEGA | 214 | ROBUST STATISTICAL METHODS | 212-226 | (7) ) STATISTIČNI h | DEFINITION OF APPROPRIATE GEODETIC GEODETSKI VESTNIK |60/2| ts = (HTEH)-1HTEui. (8) Matrika transformacije S in matrika H(HTEH)-1HTE sta ortogonalna projektorja, ki projicirata vektor opazovanj u. na dve ortogonalni komponenti: ts in u.. Če bi reševali model (7), bi iskali rešitev za ts pod pogoji metode najmanjših kvadratov, tokrat minimuma vsote kvadratov elementov vektorja u. s predhodno nastavljeno matriko uteži E. Model transformacije S podaja torej rešitev, ki minimizira evklidsko normo vektorja koordinatnih komponent točk mreže. Matrika H je datumska matrika notranjih vezi, ki vključuje vse točke geodetske mreže, ali drugače: ma- ^ trika, ki zagotavlja, da vse točke geodetske mreže definirajo datum (primer: prosta mreža). To pomeni, da je matrika E kot matrika uteži tista, ki določa, katerim točkam dodati večji oziroma manjši vpliv pri definiranju datuma. Želimo torej nastaviti tako matriko uteži E, ki bo prek matrike S definirala optimalni geodetski datum oziroma na koncu zagotavljala optimalno rešitev za vektor u.. V matematičnem smislu je primerno definiran geodetski datum pogoj za uspešno reševanje modela (7), ki je parametrični model (tudi vsak model izravnave je parametrični model). Parametrični modeli imajo v splošnem nekaj prednosti in slabosti (Hampel et al., 1986): — pojav grobih pogreškov vpliva na računane parametre, — model je samo idealizirana aproksimacija realnosti, ki pa je seveda tako matematično veliko lažje obvladljiva, SI | — parametričnemu modelu pripada dokaj enostaven stohastični model, — obstaja možnost, da je parametrični model napačno izbran, — parametrični modeli ločujejo informacije iz opazovanj na strogo strukturne in strogo slučajne spremenljivke. 3.1 Metoda največjega verjetja Pri definiranju geodetskega datuma v funkcionalnem modelu (7) se lahko pojavi problem točk, ki so se premaknile. Te točke so v modelu transformacije S neke vrste grobo pogrešena opazovanja (enačba (2)). Ta problem lahko rešimo z robustnimi statističnimi metodami, ki v splošnem minimizirajo vpliv grobo pogrešenih opazovanj v modelu. Optimalni geodetski datum lahko v tem primeru določimo z uporabo metode največjega verjetja (angl. MLE—Maximum Likelihood Estimation; Hampel et al., 1986) ali ocene M, ki omogoča določitev najverjetnejše vrednosti za iskani parameter in je tudi najpogosteje uporabljena ocena v linearnih (lineariziranih) modelih . V splošnem uporaba ocene M za rešitev obravnavanega funkcionalnega modela zahteva določitev minimuma funkcije odstopanj (Jäger et al., 2005; Berné Valero in Baselga Moreno, 2005) ali popravkov u. V obravnavanem primeru transformacije S iščemo minimum funkcije r (Jäger et al., 2005; Chen et al., 1990; Hampel et al., 1986; Rousseeuw in Leroy, 1987): P(ud4 ) = min. (9) Za ekstrem (minimum) funkcije p(ud.) velja, da je prvi odvod te funkcije po u^- enak 0. Prvi odvod definira vplivno funkcijo (Hampel et al., 1986), v tuji literaturi pogosto označeno z IF (angl. influence function). Vplivna funkcija predstavlja vpliv grobega pogreška na ocenjene vrednosti neznank v izravnavi (Hampel et al., 1986). Funkcija je oblike (Jäger et al., 2005; Berné Valero in Baselga Moreno, 2005): Aleš Marjetic, Klemen Kregar | DOLOČITEV PRIMERNEGA GEODETSKEGA DATUMA Z UPORABO ROBUSTNIH STATISTIČNIH METOD | DEFINITION OF APPROPRIATE GEODETIC DATUM USING ROBUST STATISTICAL METHODS | 212-226 | | 215 | | 60/2 | GEODETSKI VESTNIK i \ W \ PUd) n ¥(Ud ) =IF (Ugdj ) = ——- = 0- du„ (10) Ker je u^ odvisen od t5 , za določitev minimuma uporabimo pravilo posrednega odvajanja odvisne funkcije (Jäger et al., 2005) in imamo: dpuj(1)) d k (t )) v puk (15)) du) dt„ P k du) dt „ ■M")) du) dt„ = 0, (11) du) dt p kjer je: uk. — k-ti element vektorja u , k = 1...2m (št. koord. komponent), hk — k-ta vrstica matrike H, ki je dimenzije 2m x d (d— defekt datuma), p = t,, t,, a, ds — indeks, ki predstavlja transformacijske parametre med dvema koordinatnima sistemoma (geodetskima datumoma). Iz (11) izhaja pogoj ortogonalnosti (4) za določitev optimalne rešitve modela transformacije S: h H1 w(u)) v(u) ) ¥ t 2m\ = H1 u) 1 0 uj ] 0 ¥( u) ) 2 u2 0 0 2m ¥(u2 ) = H1 Eu= 0- (12) Za določitev minimuma funkcije p (enačbe 9—11) imamo pri uporabi ocene M več možnosti izbire različnih ocen (Chen et al., 1990; Jäger et al., 2005). Lahko izberemo več oblik funkcije p. : a) norma L_i oziroma prva norma vektorja u.: p(ui Hluil 1=min- Sp(u))=S\u)\=SK-hk 1 s\=min- Potem je vplivna funkcija teoretično definirana z (Jäger et al., 2005): (13) 1+1, p(uk )> 0 -1' p(u)) k )< 0 Aleš Marjetic, Klemen Kregar | DOLOČITEV PRIMERNEGA GEODETSKEGA | 216 | ROBUST STATISTICAL METHODS | 212-226 | (14) 1 STATISTIČNIH METOD | DEFINITION OF APPROPRIATE GEODETIC DATUM USING GEODETSKI VESTNIK |60/2| b) metoda Welsch (Jäger et al., 2005): p(uj ) =1 (a -a)2 - e))2 j = min., * = 2.985 Vplivna funkcija ima obliko: —(j(a ■ ak))2 W (15) (16) V literaturi (Jäger et al., 2005) je metoda izpeljana za standardizirani »popravek«, tu pa govorimo o dejanski vrednosti »popravka« (koordinatne spremembe med dvema epohama), zato predpostavljen faktor a pomnožimo s pripadajočo standardno deviacijo spremembe. Naslednja možnost, ki je v geodeziji precej razširjena predvsem pri iskanju grobih pogreškov v opazovanjih in pomeni modifikacijo metode Welsch, je: c) danska metoda (Jäger et al., 2005; Grigillo in Stopar, 2003): Pri danski metodi napišimo samo vplivno funkcijo: V(u* ) = i Eit—1 k k , __ kk -Uj , \Uj\ < c-ak k/ (c-a (17) kjer je: c — konstanta v intervalu [2, 3]. Za različne metode (a, b, c) in z upoštevanjem izraza (12) ter izraza za vplivno funkcijo y/{uk. ) posamezne metode a, b ali c imamo lahko različne načine definiranja diagonalnih elementov matrike E. Za določitev optimalne rešitve oziroma optimalne definicije geodetskega datuma uporabimo iterativni postopek prilagajanja »uteži« v diagonalni matriki E, ki jo na začetku izberemo kot enotsko matriko. E"=1 = I za it =1, (18) 2mx2m oziroma za nadaljnje iteracije Metoda Matrika E v iteraciji it > 2 L1 E=dif (VI u j—1 k} Welsch Danska E = diag i e it k I (19) E= it i it —1,k . „.it—1 1 \Uj