CORSO DI TOPOGRAFIA A A.A. 2006-2007
APPUNTI LEZIONI GPS PARTE II
Documento didattico ad uso interno - www.rilevamento.it
1
10 Inserimento di un rilievo GPS in cartografia. Problemi di trasformazione di coordinate L’uso del GPS per applicazioni geodetiche e topografiche, impone la gestione della trasformazione delle coordinate. Per inserire infatti un rilievo GPS in cartografia le coordinate WGS84 dei punti devono essere trasformate. In questo capitolo si affrontano i problemi di trasformazione di coordinate tipicamente utilizzate in ambito GPS.
10.1 Trasformazioni da coordinate geografiche a coordinate cartesiane e viceversa La relazione che lega le coordinate geografiche (latitudine, longitudine e quota ellissoidica) alle coordinate cartesiane geocentriche è dettata da semplici considerazioni trigonometriche. Nella formulazione per la trasformazione delle coordinate, si dovrà tener presente che la superficie di riferimento non è piana, ma è l’ellissoide. Troveremo quindi nelle formule i parametri che ne descrivono la geometria.
Figura10.1-1: Coordinate cartesiane e geografiche
X = (N + h ) ⋅ cosϕ ⋅ cosλ
Y = (N + h ) ⋅ cosϕ ⋅ senλ
[ (
) ]
(10.1-1)
Z = N ⋅ 1 − e + h ⋅ senϕ 2
con: ϕ,λ = latitudine e longitudine ellissoidiche del punto h = quota ellissoidica del punto
Documento didattico ad uso interno - www.rilevamento.it
2
N=
α
gran normale 1 − e 2 ⋅ sen 2ϕ e2 = quadrato dell’eccentricità dell’ellissoide Non è immediata la trasformazione inversa infatti, ricavando latitudine, longitudine e quota ellissoidica, si osserva che la formulazione della latitudine è implicita. Si risolve quindi con metodi numerici o applicando procedimenti di calcolo iterativi. Con riferimento alla 10.1-1 è immediato ricavare la longitudine dal rapporto Y/X. Y senλ = X senϕ
λ = arctan
Y X
Ricavando h dalla prima e dalla terza 10.1-1 si ottiene: h=
Z X - N ⋅ (1 − e 2 ) = −N senϕ cosϕ ⋅ cosλ
Moltiplicando i due termini dell’equazione per senϕ:
tan ϕ ⋅
X − N ⋅ senϕ = Z − (1 − e 2 ) ⋅ N ⋅ senϕ cosλ
Figura 10.1-2: Distanza dall’asse polare
Osservando che
R=
X = X2 + Y2 cosλ
Documento didattico ad uso interno - www.rilevamento.it
3
si ottiene:
ϕ = arctan
Z + N ⋅ e 2 ⋅ senϕ X2 + Y2
La formulazione della latitudine è dunque implicita e quindi può essere calcolata numericamente per tentativi oppure utilizzando uno dei numerosi metodi di calcolo numerico a convergenza che si trovano in letteratura. Tra questi, uno dei più noti per rapidità di convergenza e per precisioni ottenibili è il metodo proposto da Bencini che si basa sul calcolo di una quantità geometrica detta latitudine ridotta. R = X 2 + Y 2 (distanza dell'asse polare)
ϑ0 = arctan Z
Z R ⋅ 1 − e2
(valore di prima approssimazione della latitudine ridotta)
R ⋅ tanϑ0 a δϑ= a (correzione da apportare al valore ϑ 0 ) R 2 2 ⋅ (1 + tan ϑ0 ) − e ⋅ cosϑ0 a ⋅ 1 − e 2 + e 2 ⋅ senϑ0 −
ϑ = ϑ 0 + δ ϑ (valore corretto di seconda approssimazione della latitudine ridotta) Il procedimento termina quando il valore di ϑ si stabilizza e tende a zero; la convergenza è in genere rapidissima e già dopo la prima iterazione l'approssimazione raggiunta è sufficiente per la maggior parte della applicazioni. Determinato ϑ si calcola la latitudine con la seguente formulazione:
ϕ = arctan
tanϑ 1 − e2
Calcolate latitudine e longitudine è immediato ricavare h da una qualsiasi delle 10.1-1 h=
(
Z − N ⋅ 1− e2 senϕ
)
10.2 Richiami al concetto di Datum
L’operazione comunemente chiamata “messa in bolla” di uno strumento topografico, porta l’asse principale dello strumento stesso in posizione verticale. Tale asse diventa quindi parallelo alla tangente alla linea di forza del campo gravitazionale terrestre in quel punto ovvero perpendicolare alla superficie equipotenziale del geoide. Il geoide è dunque la superficie di riferimento rispetto alla quale il topografo esegue misurazioni, ma tuttavia, come ben noto, si tratta di una superficie di forma complessa difficilmente utilizzabile per calcoli geometrici semplici e per immediate determinazione di angoli e distanze. Risulta dunque più immediato l’utilizzo di forme geometriche semplificate, che “sostituiscono” la superficie del geoide e nel contempo permettono al topografo di utilizzare teoremi di trigonometria piana o sferica. Queste
Documento didattico ad uso interno - www.rilevamento.it
4
superfici sono l’ellissoide, la sfera e il piano. Sostituire localmente il geoide con una di queste forme, significa introdurre delle approssimazioni nei calcoli di angoli e distanze, tollerabili entro certi limiti di applicazione. Osservando la figura 10.2-1 si noti come tra la perpendicolare al geoide (V) e la perpendicolare all’ellissoide (n - scelto come superficie di riferimento) ci sia un angolo chiamato deviazione della verticale. Approssimare la superficie del geoide ad un ellissoide (sfera o piano) porta a errori di posizionamento sia in planimetria che in quota.
Figura 10.2-1: Deviazione dalla verticale
In planimetria gli scostamenti in coordinate valgono:
ε m = ϕa − ϕe
in latitudine
ε v = (λ a − λe ) ⋅ cos ϕ a
in longitudine
In altimetria, la “distanza” tra la quota misurata rispetto all’ellissoide e la quota misurata rispetto al geoide è chiamata ondulazione del geoide ( N ) ed è calcolabile come N = H-h
Documento didattico ad uso interno - www.rilevamento.it
5
Figura 10.2-2: Angolo di deviazione dalla verticale Definire un sistema di riferimento geodetico (Datum), significa scegliere la giusta geometria dell’ellissoide e il giusto orientamento affinché il valore di deviazione della verticale sia minimo e quindi siano minimi gli effetti dell’errore che si compie nel considerare la perpendicolare all’ellissoide diretta come la perpendicolare al geoide nel punto. Di conseguenza, vengono ridotti al minimo i valori di ondulazione del geoide e quindi si può considerare il valore della quota ellissoidica uguale al valore di quota ortometrica sapendo che l’errore che si compie in tale approssimazione è inferiore all’errore di misure geodetiche. Per l’Italia ad esempio è stato scelto come ellissoide di riferimento l’Internazionale di Hayford che ha le seguenti caratteristiche geometriche: α = 6378388 m e2 = 0,006722670022
Figura 10.2-3: Angolo di orientamento Monte Mario (P) Monte Soratte (Q)
Documento didattico ad uso interno - www.rilevamento.it
6
Questo ellissoide è stato orientato nel punto di emanazione di Monte Mario le cui coordinate sono state determinate mediante osservazioni a stelle fisse. In questo punto (baricentrico per la penisola italiana) l’ondulazione del geoide è nulla (N = 0). Ciò significa che in quel punto h = H e l’angolo di deviazione dalla verticale è nullo. L’ellissoide è quindi tangente al geoide in quel punto, ma può ruotare. E’ quindi indispensabile vincolarne le rotazioni imponendo e fissando un valore di Azimut (angolo tra il meridiano passante per Monte Mario) e l’arco di geodetica che congiunge Monte Mario a Monte Soratte. Definito il sistema di riferimento, è ora necessario fare in modo che gli utenti possano “collegare” i propri rilievi a tale sistema. Per questo motivo sono state predisposte le reti geodetiche, ovvero apparati costituiti da punti (materializzati in maniera ben visibile sul territorio), le cui coordinate sono note nel sistema di riferimento geodetico.
Figura 10.2-4: Rete geodetica Roma40 Riassumendo, la definizione di un sistema di riferimento geodetico è legata a tre elementi: 1. la scelta di un ellissoide 2. l’orientamento dell’ellissoide scelto 3. la definizione di una rete geodetica di appoggio
Documento didattico ad uso interno - www.rilevamento.it
7
Per l’Italia è stato scelto l’ellissoide di Hayford orientato a Monte Mario. La rete di appoggio è schematizzata in figura 10.2-4. La definizione di una rete di appoggio geodetica, permette dunque agli utenti di disporre di punti di coordinate note in un determinato Datum, sul territorio. Per inserire quindi un rilievo topografico in un Datum geodetico è sufficiente collegarsi ai vertici di tale rete. Per l’Italia esistono due sistemi geodetici fondamentali: • •
il Datum nazionale Roma40 – rete di appoggio nata per tutte le attività geodetiche, cartografiche e topografiche - ottenuta con misurazioni classiche il Datum WGS84-ETRF89 – rete di impianto IGM95, ottenuta con misure GPS.
Le due reti sono strettamente collegate tra di loro e la presenza di un numero considerevole di punti doppi ha permesso la determinazione dei valori locali dei 7 parametri, pubblicati per ciascun punto e applicabili fino a distanze di 15-20 km dal punto stesso.
Figura 10.2-5: Rete geodetica IGM95 Si osservi inoltre, che il WGS84 utilizza un ellissoide differente da quello di Hayford (GRS80) con le seguenti caratteristiche geometriche.
Documento didattico ad uso interno - www.rilevamento.it
8
α = 6378137 m e2 = 0,006694379990
Tale ellissoide non è orientato in l’Italia in quanto approssima mediamente la forma della terra a livello mondiale. Non è quindi trascurabile il valore di deviazione della verticale rispetto a questo ellissoide. Dunque diventa importante conoscere il valore dell’ondulazione del geoide N rispetto all’ellissoide WGS84 per poter ottenere dalla quota ellissoica determinata con il GPS, la quota sul livello del mare.
Figura 10.2-6: Schemi di ellissoide locale (Hayford orientato a M.M.) e ellissoide geocentrico (WGS84) Nei paragrafi successivi vedremo in dettaglio le problematiche di trasformazione tra Datum differenti e le problematiche di determinazione della quota ortometrica con misure GPS.
10.3 Trasformazione a 7 parametri per il cambiamento di Datum
Tra due sistemi di riferimento cartesiani tridimensionali, le trasformazioni di coordinate sono geometricamente governate da 6 parametri (tre traslazioni e tre rotazioni). A questi si aggiungerà un fattore di scala per le trasformazioni specifiche di Datum dal sistema WGS84 al sistema Roma40. Con riferimento alla figura 10.3-1, note le coordinate di un generico punto P nel sistema di riferimento cartesiano O1 , X1 ,Y1 , Z1 si possono determinare le coordinate del medesimo punto nel sistema di riferimento cartesiano O 2 , X 2 ,Y2 , Z 2 con una relazione del tipo: X 2 = X 0 + (1 + k ) ⋅ R ⋅ X 1
Documento didattico ad uso interno - www.rilevamento.it
(10.3-1)
9
Figura 10.3-1: Parametri di trasformazione coordinate tra due sistemi cartesiani geocentrici X1 X1 = Y1 Z1
è il vettore delle coordinate note del punto P nel sistema cartesiano O1 , X1 ,Y1 , Z1 , X2 X 2 = Y2 Z2
è il vettore delle coordinate incognite del punto P nel sistema cartesiano O 2 , X 2 ,Y2 , Z2 , X0 X 0 = Y0 Z0
Documento didattico ad uso interno - www.rilevamento.it
10
è il vettore delle componenti di traslazione ed infine cosR y cosR z
cosR x sinR z + sinR x sinR y cosR z
sinR x sinR z − cosR x sinR y cosR z
R = − cosR y sinR z
cosR x cosR z − sinR x sinR y sinR z
sinR x cosR z + cosR x sinR y sinR z
− sinR x cosR y
cosR x cosR y
sinR y
è la matrice di rotazione che si ottiene facendo ruotare la terna cartesiana con indice 1 attorno ai propri assi X1 ,Y1 , Z1 (esattamente in quest'ordine) con rotazioni pari a, rispettivamente, R x , R y , R z , per disporne gli assi parallelamente a quelli della terna cartesiana di indice 2. Vediamo ora passo per passo, il procedimento per la determinazione della matrice R analizzando separatamente le tre rotazioni e infine combinandole per ottenere il risultato complessivo. Per convenzione, si considerano di segno positivo gli angoli di rotazione in senso antiorario. Rotazione attorno all’asse X (RX)
Con riferimento alla figura 10.3-2, per trasformare le coordinate del punto P dal sistema di riferimento O, X2, Y2, al sistema di riferimento O, X1, Y1 si dovrà ruotare la terna indice 1 di un angolo pari a Rx in senso antiorario.
Figura 10.3-2: Rotazione attorno all’asse X
Y2 = AP − BC Z 2 = OC + AB
Documento didattico ad uso interno - www.rilevamento.it
11
Y2 = Y1 ⋅ cos(400 g - R X ) − Z1 ⋅ sen(400 g - R X ) = Y1 ⋅ cosR X + Z1 ⋅ senR X Z 2 = Y1 ⋅ sen (400 g - R X ) + Z1 ⋅ cos(400 g - R X ) = − Y1 ⋅ senR X + Z1 ⋅ cosR X Quindi la matrice di rotazione attorno a X:
1
0
R X = 0 cosR X 0 − senR X
0 senR X cosR X
Rotazione attorno all’asse Y (RY)
Analogamente a quanto visto nel paragrafo precedente, con riferimento alla figura 10.3-3, possiamo scrivere:
Figura 10.3-3: Rotazione attorno all’asse Y
X 2 = PA + BC Z 2 = OC − AB
Documento didattico ad uso interno - www.rilevamento.it
12
X 2 = X 1 ⋅ cos(400 g - R Y ) + Z1 ⋅ sen(400 g - R Y ) = X 1 ⋅ cosR Y − Z1 ⋅ senR Y Y2 = Z1 ⋅ cos(400 g - R Y ) − X1 ⋅ sen(400 g - R Y ) = Z1 ⋅ cosR Y + X 1 ⋅ senR Y Dalle quali è immediato ricavare la matrice di rotazione attorno a Y:
RY =
cosR Y
0 − senR Y
0 senR Y
1 0
0 cosR Y
Rotazione attorno all’asse Z (RZ) Analoghe considerazioni ci portano alla determinazione della matrice di rotazione attorno all’asse Z Con riferimento alla figura 10.3-4, possiamo scrivere:
Figura 10.3-4: Rotazione attorno a Z
Si osservi inoltre che: X 2 = OA + BC Y2 = CP − OD
Documento didattico ad uso interno - www.rilevamento.it
13
quindi X 2 = X 1 ⋅ cosR Z + Y1 ⋅ senR Z Y2 = Y1 ⋅ cosR Z − X1 ⋅ senR Z
Da queste relazioni è immediato ricavare la matrice di rotazione attorno a Z: cosR Z
R Z = − senR Z 0
senR Z
0
cosR Z 0
0 1
Il prodotto righe per colonne delle tre matrici di rotazione fornisce la matrice R complessiva:
cosR y cosR z R = R Z ⋅ R Y ⋅ R X = − cosR y sinR z sinR y
cosR x sinR z + sinR x sinR y cosR z cosR x cosR z − sinR x sinR y sinR z
sinR x sinR z − cosR x sinR y cosR z sinR x cosR z + cosR x sinR y sinR z
− sinR x cosR y
cosR x cosR y
Se le rotazioni tra i due sistemi di riferimento sono piccole e se si trascurano i termini del secondo ordine, la matrice R può essere facilmente linearizzata e restituita in forma più semplice: 1 R L = −Rz Ry
Rz
−Ry
1 −Rx
Rx 1
(10.3-2)
Il fattore di scala L’Istituto Geografico Militare ha calcolato per gli utenti, il valore dei 7 parametri di trasformazione. A ciascun vertice della rete IGM95 è associata una serie di parametri utilizzabili in un range di 15-20 km (se si vuole sfruttarne al massimo le caratteristiche di precisione: 3-5 cm). Per calcolare questi parameri, l’IGM ha utilizzato i cosiddetti “punti doppi” ovvero punti le cui coordinate sono note sia nel sistema di riferimento Roma40 che nel sistema di riferimento WGS84. Tali vertici appartengono alle due reti descritte al paragrafo precedente. Le reti, come già detto, sono state ottenute con differenti tecniche di misura e con differenti precisioni. Il fattore di scala, adatta al meglio e localmente la rete WGS84 “deformandola”sulla rete Roma40. Di fatto, questo fattore di scala modella le differenti distorsioni delle due reti geodetiche. Tenendo dunque conto del fattore di scala, la matrice R assume questa forma:
(1 + k ) Rk = − Rz Ry
Rz
− Ry
− Rx
Rx (1 + k )
(1 + k )
Documento didattico ad uso interno - www.rilevamento.it
(10.3-3)
14
10.5 Note operative per inserire un rilievo GPS in cartografia Nei paragrafi precedenti abbiamo quindi visto come la trasformazione delle coordinate ottenute nel sistema geodetico WGS84 in coordinate geografiche Roma40 sia governata da 7 parametri di trasformazione. Se alle coordinate Roma40 così ottenute si applica ora una proiezione cartografica, è immediato l’inserimento in cartografia dei punti di un rilievo GPS.
Figura 10.5-1: Schema di proiezione cartografica
Bisogna però osservare che i 7 parametri governano le trasformazioni tra il datum geodetico ETRF89 e il datum geodetico Roma40. Operativamente, ciò significa che le nostre misure geodetiche devono essere collegate alla rete IGM95 o a vertici equivalenti che possono essere passivi (vertici della rete IGM95 o di suoi raffittimenti) o vertici attivi (reti di stazioni permanenti). I primi sono materializzati sul territorio tramite chiodi topografici, i secondi dispongono invece di un GPS sempre attivo.
Figura 10.5-2: Schema di rete vincolata a due vertici WGS84 - ETRF89
Documento didattico ad uso interno - www.rilevamento.it
15
Figura 10.5-3: La rete di stazioni GPS permanenti della Regione Lombardia (attivi)
Figura 10.5-4: Vertice di raffittimento della Regione Lombardia (passivi)
Documento didattico ad uso interno - www.rilevamento.it
16
11. L'altimetria con il GPS Come gia accennato più volte i questo capitolo, il valore di quota ellissoidica rispetto all’ellissoide WGS84 ottenuto direttamente da un rilievo GPS è un parametro altimetrico puramente geometrico e spesso non collegata alla realtà fisica del campo gravitazionale terrestre. Il dislivello tra due punti, ottenuto per differenza di due quote ellissoidche non sempre è uguale al dislivello reale, ottenuto per differenza tra due quota ortometriche. In alcuni casi può essere necessario eseguire un confronto tra dati ottenuti da livellazioni geometriche, riferiti per natura alla superficie geoidica, e dati derivanti da misure GPS, invece strettamente legati all'ellissoide. Da queste considerazioni si capisce l'importanza di determinare i legami tra la quote ellissoidiche, misurate lungo la normale ad un ellissoide, e le quote ortometriche misurate lungo la verticale in un punto.
Figura 11.1: Quota ortometrica e quota ellissoidica
Con riferimento alla figura 11.1, conoscendo l'ondulazione del geoide rispetto all'ellissoide geocentrico, è possibile determinare con buona approssimazione le quote ortometriche nel seguente modo: h=H+N
(11-1)
dove: h = quota ellissoidica H = quota ortometrica N = ondulazione del geoide rispetto all'ellissoide locale H = h−N
(11-2)
Si tratta di un calcolo approssimato poiché la tangente alla linea di forza del campo di gravità della terra in un generico punto della superficie terrestre non è parallela alla normale all'ellissoide nel medesimo punto, ma, come già più volte detto, essa è inclinata rispetto a quest'ultima, di un angolo ε detto "deviazione della verticale". Per passare da un dislivello ellissoidico ad un dislivello ortometrico è necessario quindi conoscere l'ondulazione del geoide rispetto all'ellissoide considerato nei due punti P e Q.
Documento didattico ad uso interno - www.rilevamento.it
17
Le precisioni caratterizzanti queste grandezze devono essere tali per cui, il passaggio da quote ellissoidiche a quote ortometriche (o viceversa) non introduca incertezze troppo elevate per il tipo di applicazione topografica che si sta considerando.
11.1 Il software ufficiale dell’Isituto Geografico Militare per la trasformazione delle coordinate – Verto Noto il valore puntuale di ondulazione del geoide, è immediato determinare la quota sul livello del mare di un punto rispetto. E’ possibile costruire modelli matematici locali di ondulazione del geoide, purché vengano stazionati con GPS un numero sufficiente di punti di quota sul livello del mare nota. Nella realtà operativa, e per la maggior parte delle applicazioni cartografiche e topografiche, è ormai di uso comune il software ufficiale dell’IGM, Verto.
Figura 11.1-1: Quota ortometrica e quota ellissoidica
Dalle coordinate WGS84, Verto restituisce le coordinate geografiche e cartografiche nei sistemi di riferimento più utilizzati per la cartografia italiana. La quota sul livello del mare è ottenuta mediante il modello di ondulazione del geoide messo a punto dal Politecnico di Milano.
Documento didattico ad uso interno - www.rilevamento.it
18
