Topografia e cartografia digitalearch.areaopen.progettotrio.it › html › mod_prof › cartografia...

Regione Toscana

Università degli studi di Firenze

Prof. Fausto Sacerdote

Topografia

e cartografia digitale

Capitolo 4.1

Compensazione con il metodo

dei minimi quadrati

dispense del corso

Modulo Professionalizzante

Corso per Tecnico in Cartografia Tematica

per i Sistemi Informativi Territoriali

COMPENSAZIONE CON IL METODO DEI MINIMI QUADRATI

1. Introduzione

Le misure geodetiche e topografiche, che in molti casi non riguardano solo distanze e angoli, ma anche quantitanon puramente geometriche, come ad esempio l’intensita del campo della gravita, sono legate alle quantitada stimare, che in generale sono coordinate o quote di punti, da sistemi di equazioni, che costituiscono ilmodello della rete. Il modello e spesso puramente geometrico, ma puo anche contenere aspetti dinamici,come nel caso della determinazione di differenze di quote ortometriche o di geopotenziale, che richiedonomisure del modulo della gravita.

Tutte le misure sono affette da errori casuali, descritti dalla matrice di covarianza delle osservazioni, checostituisce il modello stocastico. In molti casi le osservazioni sono stocasticamente indipendenti, nel sensoche gli errori casuali di ciascuna misura non sono influenzati dagli errori casuali delle altre, e sono quindiincorrelati. La matrice di covarianza e quindi diagonale, e i suoi elementi sono le varianze delle singoleosservazioni, che indicano il livello di precisione con cui sono state eseguite. Si ricorda che la radice quadratadella varianza e detto scarto quadratico medio (sqm); spesso l’incertezza di una misura e rappresentata daun intervallo di semiampiezza uguale a 3 volte lo sqm. Questo perche nel caso di distribuzione gaussiana,che viene generalmente assunta per gli errori casuali, la probabilita che lo scarto rientri in questo intervalloe circa 0.9975, ossia molto vicina alla certezza.

Gli errori casuali di misura si trasmettono alle quantita da stimare in maniera non semplice nel caso dimodelli non lineari. Nei modelli lineari (o linearizzati, nel caso che l’approssimazione lineare sia accettabile)si puo invece applicare la regola di propagazione della covarianza: se il modello ha la forma z = Au + b(prescindendo per il momento dal significato delle grandezze z e u ), si ottiene CZ = ACUAT . Si notiche questa regola e valida indipendentemente dalle dimensioni dei vettori z e u .

In generale nelle reti geodetiche e topografiche viene eseguito un numero ridondante di misure, e quindii sistemi di equazioni di osservazione, che esprimono le quantita misurate in funzione delle grandezze dastimare, contengono piu equazioni che incognite, e in generale non possono essere risolte esattamente proprioa causa degli errori di misura, per effetto dei quali misure ridondanti possono non essere compatibili con ilmodello.

ESEMPIO: 3 dislivelli misurati a partire da una quota Q0 fissata

y1 = q1 − Q0

y2 = q2 − Q0

y3 = q2 − q1

(1)

ovvero, in forma matriciale,

⎛⎝ y1

y2

y3

⎞⎠ =

⎛⎝ 1 00 1−1 1

⎞⎠ (q1

q2

)+

⎛⎝−Q0

−Q0

0

⎞⎠ (2)

Le quote q1 e q2 sono incognite. Il sistema e ridondante: 3 equazioni (osservazioni) e 2 incognite.

NOTA: non e tuttavia possibile determinare anche Q0 . Infatti il sistema

⎛⎝ y1

y2

y3

⎞⎠ =

⎛⎝ 1 0 −10 1 −1−1 1 0

⎞⎠⎛⎝ q1

q2

Q0

⎞⎠ (3)

1

Prof. Fausto Sacerdote - Topografia e cartografia digitaleCapitolo 4.1 Compensazione con il metodo dei minimi quadrati

1

e singolare (la seconda riga e somma della prima e della terza), in accordo con il fatto che l’aggiunta di unastessa quantita a tutte le quote non cambia i dislivelli, e quindi la soluzione non e unica.

Per qualsiasi scelta di q1 e q2 in (2) e verificata l’equazione

y1 − y2 + y3 = 0 (4)

.

Di conseguenza, se le osservazioni y1 , y2 , y3 non verificano l’equazione - cosa che in generale accade, poichele misure sono affette da errore - non esiste nessuna soluzione in senso proprio q1 , q2 .

La via di uscita in questo caso consiste nella ricerca di una soluzione approssimata q1 , q2 , ad esempioadottando il criterio di rendere minima la funzione

(y1 − y1)2 + (y2 − y2)2 + (y3 − y3)2 =

=(y1 − (q1 − Q0))2 + (y2 − (q2 − Q0))2 + (y3 − (q2 − q1))2(5)

Da un punto di vista geometrico, l’equazione (2) rappresenta, al variare di q1 , q2 , un piano nello spazio3-dim y1 y2 y3 . Il punto osservato, (y1, y2, y3) , a causa degli errori di misura, non appartiene a tale piano,e il punto del piano determinato in corrispondenza di q1 , q2 e il punto di minima distanza da (y1, y2, y3) .Per evidenziare la presenza di errori di misura, spesso le equazioni di osservazione vengono scritte nella forma

y1 = q1 − Q0 + ε1

y2 = q2 − Q0 + ε2

y3 = q1 − q2 + ε3

(6)

In alternativa, e anche possibile, senza introdurre i parametri q1 , q2 , procedere direttamente alla mini-mizzazione della funzione (5) rispetto a y1 , y2 , y3 sotto la condizione y1 − y2 + y3 = 0 (vedi eq. (4)),utilizzando il metodo dei moltiplicatori di Lagrange. In questo caso tuttavia e possibile soltanto ottenerevalori corretti dei dislivelli misurati, compatibili con il modello; le quote, non essendo state introdotte nelleequazioni, non possono evidentemente essere determinate in questo stadio della elaborazione.

Il criterio adottato e simmetrico nelle tre variabili y1 , y2 , y3 ; eseguendo i calcoli risulta l’uguaglianza degliscarti |y1 − y1| , |y2 − y2| , |y3 − y3| . Questo risultato non e accettabile se e noto a priori che le tre misuresono state eseguite con differente precisione, per cui ci si aspetta scarti piu piccoli in corrispondenza dellemisure piu precise.

Una possibile strada per tenere conto di questo fatto consiste nel modificare il criterio di minimizzazioneintroducendo nella (5) dei pesi, cioe dei coefficienti moltiplicativi alle diverse componenti, tanto piu grandiquanto piu si assume che le misure delle corrispondenti componenti siano precise:

c1(y1 − y1)2 + c2(y2 − y2)2 + c3(y3 − y3)2 (7)

dove, ad esempio, si puo scegliere ci = 1/σ2i .

Generalizzando l’esempio sopra illustrato, si assume di avere a disposizione misure di diverse quantita legatefra di loro da un modello. Un altro esempio tipico in topografia sono le misure di angoli e distanze in unarete, legati fra di loro da relazioni geometriche (ad esempio un triangolo e completamente determinato dallelunghezze dei suoi lati, e quindi gli angoli sono esprimibili in funzione dei lati).

NOTA: Il metodo dei minimi quadrati puo essere applicato anche in problemi di altra natura: un caso tipico equello di misure eseguite in istanti diversi di una quantita che evolve nel tempo secondo una ben determinatalegge (ad esempio polinomiale), con coefficienti da determinare:

2


2

yi = f(ti; a0, · · · , ak) =k∑

j=0

ajtij

dove yi sono le quantita misurate (le cui misure verificano solo approssimativamente le equazioni percheaffette da errori), ti sono istanti di tempo esattamente noti, ak sono i parametri da stimare, in numeroinferiore alle quantita misurate.

A causa degli errori di misura ci si aspetta che le quantita misurate non soddisfino esattamente il modello.Ad esempio, se si misurano i 3 angoli di un triangolo, la somma delle 3 misure non e esattamente 180◦. Cisi aspetta pero che gli errori siano a media nulla, e quindi che le misure siano variabili aleatorie le cui medieverificano il modello. Il risultato che si vuole ottenere e una stima di queste medie, che vengono assuntecome valori corretti delle quantita misurate (compensazione - inglese adjustment).

Per la formalizzazione del problema, l’insieme delle n grandezze misurate (y1, · · · , yn) e descritto comeun vettore in Rn . A causa delle relazioni presenti, i valori ammissibili appartengono ad un sottoinsieme(varieta) di Rn . Le relazioni possono essere espresse da equazioni della forma fj(y1, · · · , yn) = 0 (equazionidi condizione), in numero inferiore a quello delle grandezze interessate. Ad esempio, per gli angoli αi diun triangolo,

∑i αi = 180◦ ; quindi l’insieme delle terne di angoli ammissibili e rappresentato da un piano

in uno spazio 3-dim. Se il numero delle equazioni e n − k , in realta solo k grandezze sono indipendenti; sidice allora che la varieta e k-dim.

E pero anche possibile rappresentare le relazioni fra le grandezze misurate esprimendo tali grandezze infunzione di un certo numero di parametri numerici, yi = φi(t1, · · · , tk) (nel caso delle reti topografiche sonoin generale le coordinate dei vertici della rete, eventualmente anche rotazioni, traslazioni e variazioni di scaladei sistemi di riferimento; nel caso di una grandezza che evolve sono i coefficienti del polinomio o di altraforma funzionale); se i parametri sono k , si dice allora che la varieta dei valori ammissibili e k-dim. Per unavarieta descritta da un insieme di parametri, il primo problema e proprio quello di stimare i valori numericidi questi parametri, o correzioni a valori approssimati, utilizzando le misure a disposizione; la stima dellequantita compensate viene poi fatta utilizzando le equazioni che le rappresentano in funzione dei parametri.

Puo anche accadere che siano presenti contemporaneamente parametri ed equazioni di condizione:fj(y1, · · · , yn; t1, · · · , tk) = 0 .

ESEMPI DI VARIETA

- retta nel piano y1y2:

ay1 + by2 + c = 0 oppure{

y1 = c1t + d1

y2 = c2t + d2

- circonferenza nel piano y1y2:

(y1 − a1)2 + (y2 − a2)2 = R2 oppure{

y1 = a1 + R cos αy2 = a2 + R sin α

- piano nello spazio y1y2y3:

ay1 + by2 + cy3 + d = 0 oppure

{y1 = c1u + d1v + e1

y2 = c2u + d2v + e2

y3 = c3u + d3v + e3

- sfera nello spazio y1y2y3:

3


3

(y1 − a1)2 + (y2 − a2)2 + (y3 − a3)2 = R2 oppure

{y1 = a1 + R cos α cos βy2 = a2 + R cos α sin βy3 = a3 + R sin α

- retta nello spazio y1y2y3 (intersezione di 2 piani):

{a1y1 + a2y2 + a3y3 + a0 = 0b1y1 + b2y2 + b3y3 + b0 = 0 oppure

{y1 = c1t + d1

y2 = c2t + d2

y3 = c3t + d3

- circonferenza nello spazio y1y2y3 (intersezione di 2 sfere):

{(y1 − a1)2 + (y2 − a2)2 + (y3 − a3)2 = R2

a

(y1 − b1)2 + (y2 − b2)2 + (y3 − b3)2 = R2b

con (b1 − a1)2 + (b2 − a2)2 + (b3 − a3)2 < (Ra + Rb)2 .

2. Risoluzione di problemi lineari

2.1. Stima dei parametri e delle osservabili compensate - Il problema puo ora essere espresso nellasua forma piu generale. Si assume inizialmente che le n quantita osservate, y1, · · · , yn siano fra di lorolegate in modo lineare; questo fatto viene espresso imponendo che il loro vettore rappresentativo, y ∈ Rn ,appartenga ad una varieta lineare (varieta modello) espressa introducendo sia parametri sia equazioni dicondizione:

By = At + b (8)

Se il numero di equazioni di condizione (ovvero il numero di righe della matrice B ) e m e il numerodi parametri e k (cioe t ∈ Rk ), allora la dimensione della varieta e n − m + k (deve quindi esserek < m ≤ n ). Nel caso puramente parametrico B = I , m = n ; nel caso in cui sono presenti solo equazionidi condizione A = 0 . Si assume che il rango di B e di A (quando non e nulla) sia massimo; in questocaso sia le equazioni di condizione sia i parametri sono tutti linearmente indipendenti.

A causa degli errori di misura le osservazioni yi non verificano esattamente la (8), ma possono essereespresse nella forma

By = At + b + ε (9)

Si assume invece che il vettore y sia estrazione di una variabile aleatoria Y la cui media EY verifica la(8) per un particolare valore t del parametro t .

Si tratta di determinare in funzione delle osservazioni yi una stima di tale media, ovvero un punto y ,appartenente alla varieta e quindi tale da verificare la (8) per un opportuno t , che approssimi il megliopossibile y secondo un opportuno criterio. Generalizzando i criteri di minimizzazione visti nell’esempio (cf.(5) e (6)), come funzione da minimizzare viene scelta

(y − y)T CY−1(y − y) (10)

dove CY e la matrice di covarianza di y , supposta nota a priori. Il minimo della (10) esiste certamente,mentre non esiste il massimo, dato che la (10) e continua, positiva e non limitata superiormente. Si osserviinoltre che, se la matrice CY e moltiplicata per una costante di proporzionalita, il punto di minimo della (10)

4


4

resta invariato; si puo quindi assumere che CY sia nota a meno di una costante moltiplicativa: CY = σ20Q ,

con σ20 incognito. Questa situazione corrisponde al caso realistico in cui sono note le precisioni relative,

ma non quelle assolute, delle diverse quantita misurate (ad esempio, si sa che l’accuratezza di un profilo dilivellazione e proporzionale alla radice quadrata della sua lunghezza, ma si puo non sapere la precisione diuna singola battuta).

Tenuto conto che deve essere soddisfatta la condizione By = At+b , applicando il metodo dei moltiplicatoridi Lagrange si e ricondotti alla ricerca del minimo della funzione

12(y − y)T Q−1(y − y) − (By − At − b)T λ (11)

(dove λ ∈ Rm ; il fattore 1/2 e puramente convenzionale) rispetto alle variabili y , t (quest’ultimanon compare se A = 0 , mentre, se B = I , si puo semplicemente sostituire y = At + b) nella (10) eminimizzarla rispetto alla variabile t .

Soluzione:

posto K = BQBT , N = AT K−1A , si ottiene

t = N−1AT K−1(By − b)

y = y − QBT K−1(By − b − At)(12)

In particolare:

- nel caso puramente parametrico ( B = I ):

t = (AT Q−1A)−1AT Q−1(y − b)

y = At + b(13)

- nel caso di pure equazioni di condizione ( A = 0 ):

y = y − QBT K−1(By − b) (14)

Si osservi che le stime ottenute sono funzioni lineari delle osservazioni. Questa e una conseguenza del fattoche la funzione obiettivo da minimizzare e quadratica e che la varieta modello e lineare.

Inoltre lo stimatore del parametro t e corretto. Infatti

Et = N−1AT K−1(BEy − b) = (AT K−1A)−1AT K−1(At + b − b) = t (15)

Si puo provare che fra tutti gli stimatori lineari corretti del parametro t in funzione delle osservazioni quellodefinito da (13) e lo stimatore di minima varianza (ottimalita del metodo dei minimi quadrati per modellilineari).

Si verifica anche che Ey = Ey .

NOTA: si verifica che il vettore degli scarti ε = y − y soddisfa la relazione εT Q−1At = 0 nel casoparametrico e εT Q−1v = 0 per ogni v tale che Bv = 0 nel caso di equazioni di condizione. Questerelazioni esprimono, nel caso Q = I (ovvero quando le osservazioni sono tutte indipendenti e con la stessaprecisione), l’ortogonalita degli scarti alla varieta modello. Inoltre, per quanto visto sopra, Eε = 0 .

5


5

2.2. Formule per le matrici di covarianza (ottenute per propagazione della covarianza).

- caso parametrico:

Ct

= σ20N

−1 (N = AT Q−1A)

Cε = σ20(Q − AN−1AT )

Cy

= σ20AN−1AT

(16)

- caso di eq. di condizione:

Cε = σ20QBT K−1BQ (K = BQBT )

Cy

= σ20(Q − QBT K−1BQ)

(17)

2.3. Stima di σ20 - Si verifica che, nel caso di distribuzione normale delle osservabili, la stima di massima

verosimiglianza di σ20 e

σ20 =

1n

εT Q−1ε

Questo stimatore non e pero corretto. Si verifica infatti che Eσ20 = n−k

n σ20 . Per ottenere uno stimatore

corretto σ20 basta scegliere

σ20 =

1n − k

εT Q−1ε (18)

NOTA: La stima (18) e legata alla distanza del vettore delle osservazioni dalla varieta (per Q = I e ladistanza euclidea). Il fattore 1/(n− k) dipende dal fatto che il vettore ε e vincolato a stare in una varietacomplementare alla varieta modello, che e k-dim (per Q = I , come si e visto, e la varieta ortogonale), edha quindi n − k componenti indipendenti.

3. Problemi non lineari

In generale, i problemi di compensazione che si presentano in topografia sono non lineari, e non rientranoquindi nello schema sopra illustrato.

ESEMPI:

1. (Fig.1) Si vogliono stimare le coordinate x1 , x2 su un piano di un punto P di cui si conoscono ledistanze yi da 4 punti Pi di coordinate note x

(i)1 , x

(i)2 . Le equazioni di osservazione sono

y2i − (x1 − x

(i)1 )2 − (x2 − x

(i)2 )2 = 0 i = 1, · · · , 4 (19)

Si ha quindi un sistema di 4 equazioni (tante quante le osservabili) con 2 parametri incogniti x1 , x2 , cherappresenta una varieta non lineare 2-dim in uno spazio 4-dim. Naturalmente, quando vengono introdotti irisultati delle misure, le equazioni non sono verificate esattamente a causa degli errori di misura.

2. (Fig.2) Si vuole determinare la distanza d di un punto P da una retta, avendo misurato le distanzed1 e d2 di P da 2 punti P1 e P2 sulla retta e gli angoli α1 e α2 fra i segmenti PP1 e PP2 e laperpendicolare per P alla retta. In questo caso le osservabili sono 4: d1 , d2 , α1 , α2 e le equazioni sono2, con 1 parametro d :

6


6

d − di cos αi = 0 i = 1, 2 (20)

La varieta rappresentata e 3-dim (4-2+1). Equivalentemente la varieta puo essere rappresentata dalla solaequazione di condizione (senza parametro)

d1 cos α1 = d2 cos α2 (21)

La (21) puo essere usata per trovare i valori compensati delle osservabili, da cui la stima del parametro puoessere ottenuta usando una delle (20).

OSSERVAZIONE: L’esempio 2 rappresenta una situazione particolare, che in pratica si incontra raramente,ed e qui introdotto soprattutto per illustrare la possibilita che si presentino contemporaneamente parametri edequazioni di condizione. Equazioni come quelle dell’esempio 1, invece, compaiono sempre in corrispondenzadi misure di distanza e, insieme con le equazioni di osservazione per gli angoli descritte in dettaglio nelcapitolo sulle reti topografiche, sono quelle che si incontrano abitualmente nella compensazione delle reti.

Gli esempi rientrano nella forma generale della varieta modello

F(t,y) = 0 (22)

dove F e un vettore m-dim, mentre t e y sono, come sopra, vettori rispettivamente k-dim e n-dim. Inlinea di principio niente impedisce di formulare il metodo dei minimi quadrati in modo analogo alla (11),con la stessa funzione obiettivo da minimizzare e la condizione lineare sostituita da una non lineare:

12(y − y)T Q−1(y − y) − F(t, y)

Tλ (23)

Tuttavia un’applicazione diretta porterebbe a equazioni non lineari che in generale hanno soluzione nonunica e non facile da determinare.

Di solito, tuttavia, le osservabili sono va fortemente concentrate intorno alla loro media, e quindi il loropunto rappresentativo non si allontana molto dalla varieta; questo giustifica l’approssimazione della varieta(23) con una varieta linearizzata:

F(t0,y0) +∂F∂t

(t0,y0)(t − t0) +∂F∂y

(t0,y0)(y − y0) = 0 (24)

Ovviamente il punto iniziale y0 deve essere in prossimita del valore osservato y , e il corrispondente valoreiniziale del parametro t0 deve essere tale che la (22) sia approssimativamente soddisfatta. Una possibilescelta consiste nell’imporre y0 = y e nel determinare t0 risolvendo rispetto a t k delle m componentidi (22). Questa scelta non e pero obbligata. In molti casi (ad esempio, quando si ristimano in base a nuovemisure le coordinate di punti noti) valori approssimati dei parametri sono noti a priori. In questo caso, se leequazioni di osservazione sono espresse in forma esplicita:

y = f(t) (25)

i valori approssimati delle osservabili si ottengono direttamente dalla (25).

Ci si riconduce quindi a risolvere un problema lineare col metodo dei minimi quadrati, dove, nelle formulerisolutive (12),

7


7

B → ∂F∂y

(t0,y0) ; A → −∂F∂t

(t0,y0) ; b → −F(t0,y0) +∂F∂y

(t0,y0)y0 +∂F∂t

(t0,y0)t0 (26)

Frequentemente, anziche t , y , vengono usate come incognite i loro incrementi δt = t− t0 , δy = y− y0 .La procedura risolutiva non ne risulta modificata, salvo che in luogo di (26) si ha b → −F(t0,y0) , mentreB e A non sono modificate.

Ovviamente la soluzione ottenuta appartiene alla varieta linearizzata anziche alla varieta modello non lineare,ed e quindi affetta da un errore di linearizzazione, che puo essere evidenziato se si ha una conoscenza a prioridel valore di σ2

0 e si scopre che σ20 e molto piu grande (ma, come si vedra nel capitolo sulla verifica di

ipotesi, questo puo essere dovuto anche ad altre cause).

In ogni caso, e possibile iterare la procedura usando come punto iniziale, anziche t0 , y0 , la stima ottenutat , y , e ripetere le iterazioni finche il valore di σ2

0 non si stabilizza.

NOTA: se la varieta modello e espressa in forma esplicita y = f(t) , e possibile, dopo aver ottenuto la stimat con la procedura linearizzata, ricavare le osservazioni compensate direttamente sulla varieta modello:y = f(t) . Cio non toglie che siano ugualmente presenti errori di linearizzazione che possono essere ridotticon iterazioni.

ESEMPIO (topografico)

Sia P un punto sulla facciata di un edificio, che viene rilevato facendo stazione in due punti noti P1 eP2 situati in una piazza piana. Viene definito un sistema di riferimento con l’origine in P1 , il semiasse xpositivo passante per P2 e il semiasse z positivo orientato verso l’alto nella direzione della verticale rilevatadallo strumento. Da ciascuno dei due punti di stazione si collima P e si misurano le distanze, gli angoliazimutali P ′P1P2 e P1P2P ′ , dove P ′ e la proiezione di P sul piano xy, (entrambi in verso orario) egli angoli zenitali delle direzioni di collimazione (fig.3). Per semplicita si suppone che gli angoli azimutalivengano misurati direttamente anziche come differenze di angoli riferiti allo zero dello strumento, la cuidirezione rispetto agli assi del sistema di riferimento e ignota.

Equazioni di osservazione:

d21 = x2 + y2 + z2

θ1 = arctany

x

ζ1 = arctan(x2 + y2)1/2

z

d22 = (x − x2)2 + y2 + z2

θ2 = arctany

x2 − x

ζ2 = arctan((x − x2)2 + y2)1/2

z

Equazioni di osservazione linearizzate:

8


8

d21 + 2d1δd1 = x2 + y2 + z2 + 2xδx + 2yδy + 2zδz

θ1 + δθ1 = arctany

x− y

x2 + y2δx +

x

x2 + y2δy

ζ1 + δζ1 = arctan(x2 + y2)1/2

z+

x

(x2 + y2)1/2

z

x2 + y2 + z2δx +

y

(x2 + y2)1/2

z

x2 + y2 + z2δy − (x2 + y2)1/2

x2 + y2 + z2δz

d22 + 2d2δd2 = (x − x2)2 + y2 + z2 + 2(x − x2)δx + 2yδy + 2zδz

θ2 + δθ2 = arctany

x2 − x+

y

(x2 − x)2 + y2δx +

x2 − x

(x2 − x)2 + y2δy

ζ2 + δζ2 = arctan((x − x2)2 + y2)1/2

z+

x − x2

((x − x2)2 + y2)1/2

z

(x − x2)2 + y2 + z2δx+

+y

((x − x2)2 + y2)1/2

z

(x − x2)2 + y2 + z2δy − ((x − x2)2 + y2)1/2

(x − x2)2 + y2 + z2δz

Come valori approssimati si possono scegliere ad esempio per le osservabili gli stessi valori misurati: di =di , θi = θi , ζi = ζi , per le incognite x , y , z i valori ottenuti dalle prime 3 equazioni:

x = d1 sin ζ1 cos θ1

y = d1 sin ζ1 sin θ1

z = d1 cos ζ1

I termini noti vengono tutti portati al secondo membro e, in virtu della scelta dei valori approssimati, devonorisultare nulli per le prime 3 equazioni, mentre gli altri 3 devono essere piccoli (di un ordine di grandezzacoerente con l’accuratezza delle misure). Il verificarsi di queste condizioni puo essere utilizzato come primocontrollo che non siano stati fatti grossolani errori di calcolo. Conviene inoltre dividere le equazioni delledistanze per 2di , per ricondurre il modello lineare generale, By = At+b alla forma puramente parametrica,B = I ; in caso contrario B e una matrice quadrata diagonale, ma gli elementi corrispondenti alle equazionidelle distanze non sono uguali a 1.

In conclusione si ottiene un modello puramente parametrico, con

A =

⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

xd1

yd1

zd1

− yx2+y2

xx2+y2 0

x(x2+y2)1/2

zx2+y2+z2

y(x2+y2)1/2

zx2+y2+z2 − (x2+y2)1/2

x2+y2+z2 δzx−x2

d2

yd2

zd2

y(x2−x)2+y2

x2−x(x2−x)2+y2 0

x−x2((x−x2)2+y2)1/2

z(x−x2)2+y2+z2

y((x−x2)2+y2)1/2

z(x−x2)2+y2+z2 − ((x−x2)

2+y2)1/2

(x−x2)2+y2+z2 δz

⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

b =

⎛⎜⎜⎜⎜⎜⎜⎜⎝

000

(x−x2)2+y2+z2−d2

22d2

arctan yx2−x − θ2

arctan ((x−x2)2+y2)1/2

z − ζ2

⎞⎟⎟⎟⎟⎟⎟⎟⎠Come si e gia visto, in generale la matrice C e diagonale, e i suoi elementi sono le varianze delle osservazioni.Assumendo ad esempio che lo sqm delle distanze sia di 2mm, quello degli angoli azimutali 3arcsec, quellodegli angoli zenitali 5arcsec, usando come unita di misura per le distanze il metro, per gli angoli il radiantesi ha σ2

d = 4 ∗ 10−6 , σ2θ = 2.12 ∗ 10−10 , σ2

ζ = 5.88 ∗ 10−10 . Scegliendo σ20 = 10−10 , Si ottiene

9


9

Q =

⎛⎜⎜⎜⎜⎜⎝40000 0 0 0 0 0

0 2.12 0 0 0 00 0 5.88 0 0 00 0 0 40000 0 00 0 0 0 2.12 00 0 0 0 0 5.88

⎞⎟⎟⎟⎟⎟⎠Quando si eseguono linearizzazioni, e in generale quando nelle equazioni compaiono angoli e distanze, gliangoli vanno obbligatoriamente misurati in radianti (fa eccezione il caso di problemi lineari con soli angoli,per cui l’unita di misura e arbitraria); l’unita di misura delle distanze e in ogni caso arbitraria, basta esserecoerenti nel corso dell’intero calcolo. E opportuno evitare scelte di unita di misura che portino a moltiplicarenumeri enormemente grandi per numeri enormemente piccoli, per ragioni di stabilita numerica. Non enecessaria un’elevata precisione per gli elementi della matrice Q , che possono quindi essere arrotondati.

Nella formula risolutiva, t = (AT Q−1A)−1AT Q−1(y−b) , t e y sono in questo caso incrementi rispettoai valori apporssimati; avendo scelto i valori approssimati delle osservabili uguali ai valori misurati, si deveporre y = 0 .

Lo schema qui adottato, a parte qualche semplificazione dovuta alla particolare scelta del sistema di rifer-imento, puo essere facilmente generalizzato ad una qualsiasi rete topografica, per cui quindi e abbastanzafacile tradurre questa procedura in un programma di calcolo di uso generale.

4. Verifiche di ipotesi (varieta lineare espressa parametricamente: y = At + b )

Le stime ai minimi quadrati possono essere utilizzate per verifiche di ipotesi, applicando test parametriciper accettare o rifiutare informazioni a priori sul valore di certi parametri. Le applicazioni piu significativeriguardano la precisione globale delle misure, espressa dal parametro σ2

0 , e gli stessi parametri ti (adesempio, le coordinate dei vertici di una rete), di cui possono essere noti valori ottenuti in campagne dimisura diverse da quella su cui e stata eseguita la compensazione. Particolare rilevanza hanno i test dicontrollo delle deformazioni, in cui i risultati di una campagna di misure eseguite ad esempio su un manufatto(edificio, diga, ecc.) vengono confrontati con i risultati di campagne precedenti, assumendo come ipotesi daverificare che gli spostamenti siano nulli (ovvero che le differenze dei risultati siano puramente dovute aglierrori casuali delle misure).

Per poter eseguire un test parametrico e necessario conoscere il tipo di distribuzione di probabilita dellastima utilizzata nel test. Si premettono quindi alcune considerazioni sulle distribuzioni di probabilita dellequantita che intervengono nelle stime ai minimi quadrati.

. y (n-dim) vettore delle osservabili (con matrice di covarianza C = σ20Q ). Si assume sempre che la sua

densita di probabilita congiunta sia gaussiana

. t (k-dim) vettore dei parametri

Si puo provare che, se y e gaussiana e C e la sua matrice di covarianza, (y − y)T C−1(y − y) hadistribuzione χ2

n ( y = Ey = Ey )(vedi Appendice).

Inoltre (σ20)−1(t−t)T N(t−t) = (y−y)T C−1(y−y) hanno distribuzione χ2

k (Ricordare che N = AT Q−1Ae che σ2

0N−1 e la matrice di covarianza di t : inoltre y = At + b , y = At + b) .

Di conseguenza, si puo provare che (y− y)T C−1(y− y) ha distribuzione χ2n−k . (Puo servire ricordare che

(y− y)T C−1(y− y) = (y− y)T C−1(y− y)+ (y− y)T C−1(y− y) e tener presente la regola della somma didistribuzioni χ2 indipendenti. Per verificare l’uguaglianza basta tener presente che y−y = (y−y)+(y−y) ,e che (y − y)T C−1(y − y) = 0 ; quest’ultima uguaglianza puo essere provata semplicemente sostituendo ay la sua espressione in funzione delle osservabili).

NOTA: In generale queste statistiche non possono essere usate direttamente per i test, poiche contengono il

10


10

parametro σ20 , che di solito e incognito.

4.1. Test su σ20 - Il test su σ2

0 serve per verificare l’attendibilita di un’informazione a priori sull’accuratezzaglobale delle misure.

Ricordare σ20 = (n−k)−1εT Q−1ε (ε = y−y) . Per quanto visto sopra, (n−k)(σ2

0/σ20) = (y−y)T C−1(y−y)

ha distribuzione χ2n−k . Si vuole verificare un’ipotesi H0 : σ2

0 = σ20 . Dato il livello di significativita α ,

poiche in generale si e interessati a rifiutare l’ipotesi quando il valore di σ20 e troppo grande (il vettore delle

osservazioni si discosta troppo da quello delle osservabili compensate), si prende un intervallo di accettazionedel tipo ]0, bα[ , dove bα e tale che P ((n− k)(σ2

0/σ20) ≥ bα) = P (σ2

0 ≥ bασ20/(n− k)) = α (test a 1 coda).

Se pero si vogliono scartare anche valori piccoli di σ20 , si sceglie un intervallo del tipo ]aα/2, bα/2[ , dove

P ((n − k)(σ20/σ2

0) ≤ aα/2) = P ((n − k)(σ20/σ2

0) ≥ bα/2) = α/2 .

La presenza di un σ20 piu grande di quello che ci si aspetta puo significare diverse cose:

- effetti di linearizzazione. In questo caso pero si puo iterare la procedura di stima m.q. finche il valore diσ2

0 non si stabilizza. Se e ancora troppo grande, bisogna cercare altre cause;

- errori di modello: e sbagliata la varieta lineare a cui si impone che le osservazioni compensate appartengano.Si pensi, ad esempio, ad una rete di livellazione vincolata alle quote considerate note di 2 punti, di cui perouna e sbagliata; in questo caso, nelle equazioni di osservazione compare un termine noto errato;

- errori grossolani di misura (outliers). Assumendo che l’outlier sia uno solo, si puo tentare di individuarlocon la seguente procedura:

. si determinano gli scarti standardizzati (cioe le differenze |yi − yi| divise per i loro sqm) e si elimina lamisura corrispondente a quello piu grande;

. si esegue una nuova stima m.q. usando le restanti osservazioni (purche il loro numero continui ad essereridondante) e si determina un nuovo σ2

0 , che viene utilizzato per un nuovo test

. se l’ipotesi H0 viene accettata, si procede determinando, sulla base dei valori stimati dei parametri, ilvalore compensato dell’osservazione scartata, e quindi il nuovo scarto relativo a questa osservazione, che puoa sua volta essere utilizzato per un test, essendo distribuito come una gaussiana a media nulla e con varianzadeterminabile per propagazione.

Non sempre, tuttavia, questa procedura da risultati soddisfacenti, dato che una misura errata influenza lacompensazione di tutte le misure, e puo determinare scarti grandi anche in corrispondenza di misure corrette(si pensi, come caso limite, alla misura dei 3 angoli di un triangolo: un grosso errore in uno solo degli angolimisurati si equiripartisce nella compensazione su tutti e 3 gli angoli; quindi gli scarti risultano tutti ugualie non e possibile individuare l’outlier). Pensando ad un altro esempio concreto, una rete di livellazione, cisi rende conto anche intuitivamente che e tanto piu facile individuare l’outlier quanto piu numerose sono lemisure di dislivello che interessano i punti collegati dalla misura errata (ridondanza locale).

Bisogna inoltre tener conto che, nel caso di reti con un elevato numero di quantita misurate, e molto facile chegli outliers siano piu di uno, e vanno quindi costruite procedure iterative per una loro completa rimozione.

4.2. Test sui parametri - Per verificare un’ipotesi sul valore del parametro (vettoriale) t , si tengapresente che sia (σ2

0)−1(t− t)T N(t− t) sia (n−k)σ20/σ2

0 , di cui si conosce la distribuzione (sono entrambeχ2 ), hanno entrambe il parametro incognito σ2

0 a denominatore; il loro rapporto quindi non contiene taleparametro. Si verifica facilmente (vedi Appendice) che k−1(t − t)T N(t − t)/σ2

0 ha distribuzione Fk,n−k .Anche in questo caso in generale si rifiuta l’ipotesi quando il valore empirico di questa quantita e troppoelevato e si fa quindi un test a 1 coda.

In modo analogo si puo fare il test non su tutti i k parametri, ma solo su r fra loro. In questo caso laforma quadratica a numeratore contiene solo le r componenti prese in considerazione del vettore t − te l’inversa della loro matrice di covarianza che, si badi bene, si ottiene cancellando dalla matrice N−1 lerighe e le colonne che non riguardano le componenti in esame e invertendo la matrice cosi ottenuta (non e lastessa cosa che cancellare righe e colonne della matrice N ); la distribuzione della statistica cosi costruita e

11


11

Fr,n−k .

Se poi il parametro su cui si fa il test e uno solo, si tenga conto, sempre con riferimento all’Appendice, che[(ti − ti)/

√(N−1)

ii]/σ0 ha distribuzione t(n−k) ( (N−1)ii e la varianza di ti ). In questo caso, essendo

la densita di probabilita t simmetrica rispetto all’origine, conviene prendere un intervallo di accettazionesimmetrico.

NOTA: sia t un parametro 2-dim (ad esempio la coppia di coordinate planimetriche di un punto); fissatot (cioe il valore ipotizzato), l’insieme dei valori stimati t per cui l’ipotesi viene accettata verifica unadisequazione della forma (t − t)T C−1

t(t − t)/(kσ2

0) < bα . Si consideri ora questa stessa disequazione

tenendo fissa la stima t e facendo variare t ; si ottiene cosi l’insieme dei valori ipotizzati che potrebberoessere accettati dal test in base a quella stima. Nel piano della variabile 2-dim t questo insieme e l’internodi un’ellisse con centro in t (ellisse di confidenza). E proprio questa ellisse che viene disegnata (in unascala opportuna) nella rappresentazione grafica di una rete compensata per evidenziare l’accuratezza dellestime ottenute, e quindi la bonta del progetto di rete (ellissi grosse significano che si possono commetteregrossi errori di stima, ellissi molto allungate significano che gli errori sono particolarmente grandi in unadirezione). Si noti che in questo ragionamento σ2

0 e considerata come una costante, e influisce sullagrandezza assoluta delle ellissi, ma non sulla loro forma, che esprime l’accuratezza relativa delle coordinatestimate ed e completamente determinata dalla matrice C

t; e quindi possibile determinare a priori (sulla base

del modello di rete, prima di eseguire le misure), a meno di un fattore globale di scala, le ellissi corrispondentiai diversi vertici. Ovviamente nel caso 3-dim (coordinate plano-altimetriche di punti) si hanno ellissoidi alposto delle ellissi.

4.2. Test di controllo della stabilita - La stabilita di una struttura puo essere verificata medianteuna rete con i vertici in parte appartenenti alla struttura, in parte al di fuori di essa in posizioni tali daavere ragionevole certezza della loro stabilita; vengono eseguite campagne di misura in epoche diverse, e sisottopone a test l’ipotesi che siano nulle le differenze fra le misure delle coordinate ottenute in campagnediverse. Se l’ipotesi e verificata, le differenze di coordinate ottenute sono attribuite ad errori casuali dimisura; in caso contrario vengono interpretate come dovute a deformazioni.

Si assuma che nelle due campagne vengano eseguite esattamente le stesse misure di dislivello, secondo lostesso schema di rete, in modo che le stime dei parametri ottenute nelle due campagne, t1 , t2 siano basatesu misure indipendenti delle stesse quantita. In tal caso la differenza t2 − t1 , in assenza di deformazioni,ha media nulla e matrice di covarianza 2σ2

0N−1 ; se si vogliono usare le stime di σ2

0 ottenute da entrambele serie di misure, si ha per (n − k)(σ2

01 + σ202)/σ2

0 , che e somma di due variabili χ2n−k indipendenti, una

distribuzione χ22(n−k) .

Si usa quindi come statistica per la verifica dell’ipotesi di stabilita la variabile aleatoria

12k

(t2 − t1)T N(t2 − t1)12 (σ2

01 + σ202)

.

In modo analogo a quanto visto in relazione al test sul valore dei parametri, si puo provare che, se e verificatal’ipotesi H0 : t2−t1 = 0 , tale variabile ha distribuzione Fk,2(n−k) ; quindi, fissato un livello di significativita,si puo determinare un intervallo entro il quale le differenze di coordinate possono essere attribuite ad erroricasuali e non a deformazioni.

12


12


13


14


15


16


17


18


19


20


21


22

Minimi quadrati 1

Sistema ridondante di equazioni

Esempio

y3

y1q1

q2

y2

Q0123

022

011

qqy

Qqy

Qqy

ossia

011

10

01

0

0

2

1

3

2

1

Q

Q

q

q

y

y

y

3,2,1, iyi dislivelli misurati

2,1, jq j quote incognite

0Q quota nota


23

Minimi quadrati 2

Esiste soluzione solo se

0321 yyy (equazione di condizione)

non verificata a causa degli errori di misura

Ricerca di soluzione approssimata 21ˆ,ˆ qq

tale che, posto0

ˆ

ˆ

11

10

01

ˆ

ˆ

ˆ

0

0

2

1

3

2

1

Q

Q

q

q

y

y

y

sia minimizzato2

33

2

22

2

11 )ˆ()ˆ()ˆ( yyyyyy

3

2

1

ˆ

ˆ

ˆ

y

y

y

è la proiezione ortogonale di3

2

1

y

y

y

sul piano di equazione 0321 yyy


24

Minimi quadrati 3

Metodi di risoluzione

si sostituisce)(

ˆˆˆ

ˆˆ

ˆˆ

123

022

011

qqy

Qqy

Qqy

e si minimizza rispetto a 21ˆ,ˆ qq

si applica il metodo dei moltiplicatori di Lagrange:

si minimizza

)ˆˆˆ()ˆ()ˆ()ˆ( 321

2

33

2

22

2

11 yyyyyyyyy

rispetto a 321ˆ,ˆ,ˆ yyy

sotto la condizione 0ˆˆˆ321 yyy

e si ricava da21ˆ,ˆ qq )(

soluzione:

se 321 yyy ,

risulta3

ˆ,3

ˆ,3

ˆ332211 yyyyyy

(equiripartizione dell’errore di chiusura)


25

Minimi quadrati 4

Il procedimento descritto è simmetrico rispetto agli scarti

332211ˆ,ˆ,ˆ yyyyyy

Se le misure hanno differente precisione, vanno introdotti dei pesi:

si minimizza 2

333

2

222

2

111 )ˆ()ˆ()ˆ( yycyycyyc

maggiore precisione peso maggiore

ad esempio 2

1

k

kc


26

Minimi quadrati 5

Caso generale (lineare)

bAtBy

nRy osservabili (note)

kRt parametri (incogniti)

nmB

nmk

kmA

matrice di covarianza nota, C

eventualmente a menodi un coeff. di proporzionalità

QC 2

0

incognito noto


27

Minimi quadrati 6

Si minimizza

)()ˆˆ()ˆ()ˆ(2

1 1 mTTRbtAyByyQyy

rispetto a yt ˆ,ˆ

sotto la condizione btAyB ˆˆ

Soluzione:

AKANBQBKtAbByKQByy

bByKANt 1

1

11

)ˆ(ˆ

)(ˆTT

T

T

NOTA btAyBE

ossia il valore atteso del vettore delle osservabili

sta sulla varietà modello

Inoltre ttE

ossia t è uno stimatore corretto del parametro t

e yy EE ˆ


28

Minimi quadrati 7

Casi particolari:

(equazioni di condizione0A

senza parametri)

)(ˆ 1bByKQByy

T

(caso puramente parametrico)IB

si può sostituire

btAy ˆˆ

e minimizzare

)ˆ()ˆ( 1yyQyy

T

rispetto a t

AQANbtAy

byQANt 111

ˆˆ

)(ˆT

T


29

Minimi quadrati 8

Le matrici di covarianza

di yyyt ˆ,ˆ,ˆ

si determinanocon la legge di propagazione:

se pMuv

( variabili aleatorie, qs

RuRv ,

M matrice qs ,s

Rp

con elementi numerici noti)

allora

TMMCC uv

ad esempio

12

0ˆ NCt


30

Minimi quadrati 9

NOTA le soluzioni yt ˆ,ˆ

sono indipendenti da . 2

0

La quantità )ˆ(1

ˆ 12

0 yyQT

kn

è uno stimatore corretto di2

0

(ossia ) 2

0

2

0ˆE

NOTA forma quadratica2

0ˆ

definita positiva in

2

0ˆ grande || grande

qualche (osservazione compensata)y

si discosta molto dalla misura y

errori di misura significativi


31

Minimi quadrati 10

Modello non lineare

0y)F(t,

(sistema m-dim. di di eq. di osservazione)

Linearizzazione

0yyyty

Fttyt

t

FytF ))(,())(,(),( 00000000

b A B

occorre conoscere

i valori approssimati 00 , yt

si può porre yy 0

(vettore delle osservazioni)

per ricavare il vettore k-dim 0t

(se non è noto a priori)

bisogna risolvere un sistema di k eq. scelte fra le m del sistema 0ytF ),(


32

Minimi quadrati 11

Test statistico su2

0

(entità degli errori di misura)

si assume che la distribuzione di prob.degli errori sia gaussiana

stime dei parametri gaussianet

osservazioni compensate gaussianey

scarti gaussiani

2

0

2

0ˆ

)( knX ha distribuzione2

kn :

1)2/(2/)( nx

X xCexf


33

Minimi quadrati 12

Test statistico su2

0

Si ipotizza un valore di2

0

(indice globale dell’entità degli errori)

2

0

2

0ˆ

troppo grande

errori di misura più grandi di quanto ci si aspettasse a priori

)Pr( xX

(livello di significatività)

Se è piccolo (es.: 05.0 )

è poco probabile che xX

se si trova xknX2

0

2

0ˆ

)(

con ipotizzato2

0

e ricavato dalle misure2

0ˆ

è molto probabile ( 95.01Pr )

che l’ipotesi su sia inattendibile 2

0


34

Minimi quadrati 13

Test statistico su2

0

Un valore troppo grande di2

0ˆ

può essere dovuto a:

effetti di linearizzazione

errori di modello

errori grossolani di misura (outliers)

In caso di outliers

si determinano gli scarti standardizzati (ossia divisi per i loro sqm)

si elimina l’eq. di osservazione con lo scarto standardizzato più elevato e si procede di nuovo alla stima m.q.

si verifica se in questo caso

il valore ipotizzato di2

0

risulta accettabile dal test

Si possono eliminare una per volta più equazioni di osservazione, purchè rimanga la ridondanza


35


36

Topografia e cartografia digitalearch.areaopen.progettotrio.it › html › mod_prof › cartografia...

Documents

Transcript of Topografia e cartografia digitalearch.areaopen.progettotrio.it › html › mod_prof › cartografia...