Come risolvere l’intersezione 3D da stazione totale

Gianni Rossi

Credo non ci sia un geometra topografo che non conosca l’intersezione in avanti per determinare un punto inaccessibile facendo due stazioni con la TS e misurando soltanto gli angoli al punto incognito. Il caso classico è quello di rilevare la croce sulla sommità di un campanile (spesso un trigonometrico). Mi ricordo che a scuola di geometra (almeno ai miei tempi) era uno dei primi schemi che ci insegnavano nel corso di topografia.

Al giorno d’oggi è uno schema già meno indispensabile per chi possiede la stazione totale dotata della modalità “senza-prisma”. Tuttavia, anche per chi dispone di questa tecnologia l’applicazione dell’intersezione in avanti è comunque utile per controllare la bontà della rilevazione. Questo perché, pur rilevando il punto con il senza-prisma, vengono comunque registrati gli angoli azimutale e zenitale e si è quindi in grado di applicare lo schema e ottenere così un secondo risultato che garantisca la correttezza della rilevazione diretta.

Naturalmente finora ho parlato di due sole stazioni per identificare lo schema minimo indispensabile a risolvere il problema, ma la buona prassi topografica della ridondanza di misure impone di farne almeno tre, realizzando così due intersezioni dalle quali poter avere la verifica sul risultato. Del resto l’obbligo di fare tre stazioni è imposto dalla stessa normativa catastale nel caso l’intersezione venga adottata in lavoro presentato con Pregeo, come descritto in Figura 1.

Figura 1La corretta prassi topografica prevede la doppia intersezione in avanti realizzata mediante tre stazioni, così da conseguire ridondanza di misure ed avere un controllo sul risultato. Questo vincolo è addirittura imposto dalla normativa catastale per i lavori presentati con il software Pregeo.

Dopo questa premessa vi domanderete perché nel titolo parlo di intersezione 3D e non dell’intersezione in avanti appena citata. È perché in questo articolo illustro un calcolo diverso da quello classico studiato a scuola. Quello è sostanzialmente un’algoritmo planimetrico mentre questo che sto per proporvi è effettivamente 3D. Un’altra differenza tra i due schemi è che nell’intersezione in avanti classica le due stazioni devono osservarsi reciprocamente, mentre in quella 3D questo requisito non è richiesto, a patto ovviamente che le stazioni siano determinate dalla parte di rilievo precedente, ad esempio per essere punti GPS o determinate da altre stazioni TS. Non cambiano invece le prescrizioni previste dal Catasto per realizzare l’intersezione mediante una geometria che garantisca la miglior precisione possibile. Con riferimento alla Figura 2, devono cioè essere rispettate queste condizioni:

\[ dA >= \frac{2}{3}dX\hspace{20mm}dB >= \frac{2}{3}dX \]
\[ 35^g < α < 135^g\hspace{20mm}35^g < β < 135^g \]

Figura 2Le prescrizioni previsti dalla buona tecnica topografica (e dalla normativa catastale) per una corretta intersezione a doppia base:

  • la lunghezza delle due basi misurate, dA e dB non deve essere inferiore ai 2/3 della distanza incognita dX;
  • l’angolo sul punto inaccessibile, α e β, deve essere compreso tra 35g e 165g.

Le analogie tra i due algoritmi terminano qui, per cui d’ora in avanti vedremo soltanto l’intersezione 3D. Per quelli di voi che desiderano invece approfondire l’intersezione in avanti classica, non sarà certo difficile reperire qualche documento (anche online) che la spiega[1].

Come detto, lo schema è utile quando, non disponendo di una stazione totale con misurazione senza-prisma, il punto da determinare non è raggiungibile perché interno a una proprietà recintata a cui non si ha accesso, come illustrato in Figura 3 (a sinistra), oppure nel caso classico di rilevare l’asse di un campanile collimando la croce sulla sommità (Figura 3 a destra).

Figura 3I casi di applicazione dell’intersezione 3D. A sinistra, si deve rilevare un punto all’interno di una proprietà non accessibile. A destra, si deve rilevare l’asse di un campanile collimando la croce sulla sommità.

Lo schema si chiama intersezione 3D perché permette di rilevare sia planimetricamente che altimetricamente il punto inaccessibile mediante le sole letture angolari da due stazioni, anche, come detto, senza che le stesse si osservino reciprocamente, purché siano determinabili dal resto del rilievo. Si tratta di trovare l’intersezione nello spazio delle rette 3D dei raggi di mira di ciascuna osservazione fatta dalle due stazioni, come illustrato in Figura 3. Ora, come è facile immaginare, nella realtà è praticamente impossibile che i due raggi di mira si intersechino in quanto basta anche uno scostamento impercettibile, tipo un decimo di millimetro, perché ciò non avvenga. Dobbiamo infatti pensare che nel piano, se due rette non sono parallele, si intersecano sempre, mentre nello spazio due rette non parallele possono non intersecarsi mai, anzi, questa è proprio la condizione più probabile. La Figura 4 (in alto) mostra, esasperandola, la mancata intersezione dei due raggi di mira.

Figura 4I raggi di mira delle due stazioni non si intersecano nello spazio, occorre quindi trovare l’intersezione apparente data dal punto che si trova alla minima distanza da entrambi i raggi. Tale punto e la mezzeria del segmento piu corto che connette le due rette.

Come si fa allora a trovare il punto collimato da entrambi i raggi di mira se questi non si intersecano?

Qui entra in gioco il concetto di “intersezione apparente”, vale a dire il punto che meglio di tutti approssima quello che si avrebbe se le due rette si intersecassero effettivamente. L’intersezione apparente è il punto alla minima distanza (uguale) da entrambe le rette. Come illustrato in Figura 4 (in basso), dobbiamo pertanto trovare il segmento più corto che le connette. Il punto medio di questo segmento è quello che meglio approssima il punto collimato. Sempre con riferimento alla Figura 4, è abbastanza intuitivo rilevare che il segmento minimo che congiunge le due rette è quello rosso perpendicolare ad entrambe. In Figura 4 la perpendicolarità è indicata con il classico quadratino (di colore gricio) formato dal segmento di connessione e ciascuna retta. Qualsiasi altro segmento non perpendicolare a tutte e due è palesemente più lungo. Questa doppia perpendicolarità è la proprietà che permette di risolvere il problema. Ma vediamo il tutto in dettaglio schematizzando gli elementi in gioco mediante la simbologia di Figura 5.

Figura 5A differenza che sul piano, nello spazio due rette, anche se non sono parallele, possono non intersecarsi, anzi, questa è proprio la condizione della doppia osservazione di un punto da due stazioni TS: basta anche uno scostamento minimo (frazione di mm) perché i due raggi di mira non si incrocino. Data questa condizione, il punto di intersezione si trova sul punto medio del segmento più corto che unisce le due rette, cioè quel segmento che è perpendicolare ad entrambe.

Innanzi tutto, come si evince dalla Figura 5 dobbiamo definire le due rette 3D costituite dai due raggi di mira che escono dalla TS posta nelle due stazioni. Per fare questo ci servono due punti nello spazio, il primo è ovviamente il centro strumentale della TS per entrambe le stazioni, cioè P1 per la prima stazione e P3 per la seconda. Il secondo punto di ciascuna retta invece ci manca, ma non è un problema calcolarlo. Ci basta simulare la battuta da ciascuna stazione a un punto sufficientemente distante (P2 dalla stazione P1 e P4 dalla stazione P3) in modo da calcolare in maniera precisa la retta. A questo scopo torniamo a quanto detto all’inizio su quanto previsto dalla buona tecnica topografica (e la normativa catastale) per lo schema dell’intersezione in avanti, e cioè che le due stazioni siano distanti tra loro non meno dei 2/3 della distanza al punto incognito. Sulla base di questo presupposto possiamo ragionevolmente assumere proprio la distanza tra le due stazioni quale valore per calcolare il secondo punto di ciascun raggio di mira. Troviamo pertanto la distanza tra le stazioni P1 e P3 applicando semplicemente il teorema di Pitagora alle differenze di coordinate:

\[ D_{13}=\sqrt{(x_3-x_1)^2+(y_3-y_1)^2} \]

Fatto ciò, dobbiamo determinare le coordinate 3D di questi due ipotetici punti P2 e P4 battuti idealmente dalle stazioni P1 e P3. Per farlo dobbiamo però considerare che il raggio di mira (cioè la retta su cui andremo a operare) non parte dal punto a terra ma dal centro della TS. Quindi, mentre le coordinate x e y delle stazioni P1 e P3 sono quelle derivanti dal calcolo della parte di rilievo che le ha determinate, per quanto riguarda la quota z, dobbiamo sommare l’altezza strumentale hs della TS in entrambe le stazioni:

\[ z_1= z_{1\hspace{1mm}terra}+hs_1\hspace{20mm}z_3= z_{3\hspace{1mm}terra}+hs_3 \]

A questo punto possiamo calcolare le coordinate x, y, z, dei due punti fittizi P2 e P4 facendo finta che siano stati battuti dalle rispettive stazioni P1 e P3 leggendo lo stesso angolo verticale letto al punto di intersezione da trovare (φ1 e φ3). Naturalmente per calcolare le coordinate planimetriche x e y dobbiamo conoscere l’azimut tra la stazione e tale punto di intersezione. Questo è facilmente ricavabile a partire dall’angolo orizzontale letto effettivamente dalla stazione al punto (α12 in P1 e α34 in P3) sottraendogli la correzione angolare della stazione stessa derivante dal calcolo del rilievo che l’ha determinata (corr1 in P1 e corr3 in P3):

\[ ϑ_{12}=α_{12}-corr_1\hspace{20mm}ϑ_{34}=α_{34}-corr_3 \]

Va da sé che quando si esegue la sottrazione tra due angoli, qualora il risultato fosse negativo, il valore va riportato in positivo sommando l’angolo giro (400g). Con i dati di cui sopra possiamo ora calcolare le coordinate x, y, z dei punti P2 e P4 come segue:

\[ x_2=x_1+D_{13}\hspace{1mm}sinϑ_{12}\hspace{10mm}y_2=y_1+D_{13}\hspace{1mm}cosϑ_{12}\hspace{10mm}z_2=z_1+D_{13}\hspace{1mm}cosφ_1 \]
\[ x_4=x_3+D_{13}\hspace{1mm}sinϑ_{34}\hspace{10mm}y_4=y_3+D_{13}\hspace{1mm}cosϑ_{34}\hspace{10mm}z_4=z_3+D_{13}\hspace{1mm}cosφ_3 \]

Giunti a questo punto, abbiamo a che fare con quattro punti 3D formati da due coppie che individuano ciascuna la retta nello spazio del raggio di mira emesso dalla rispettiva stazione. Da qui in avanti il processo prosegue applicando le proprietà del calcolo vettoriale. La prima di queste operazioni è quella di trovare le differenze di coordinate d13, d21, d43 tra i punti P3P1, P1P2 e P3P4:

\[ d_{13x}=x_1-x_3\hspace{10mm}d_{21x}=x_2-x_1\hspace{10mm}d_{43x}=x_4-x_3 \]
\[ d_{13y}=y_1-y_3\hspace{10mm}d_{21y}=y_2-y_1\hspace{10mm}d_{43y}=y_4-y_3 \]
\[ d_{13z}=z_1-z_3\hspace{10mm}d_{21z}=z_2-z_1\hspace{10mm}d_{43z}=z_4-z_3 \]

Poi vanno calcolati ben cinque fattori f dati dalla somma dei prodotti incrociati di tali differenze:

\[ f_{1343}=d_{13x}\hspace{2mm}d_{43x}+d_{13y}\hspace{2mm}d_{43y}+d_{13z}\hspace{2mm}d_{43z} \]
\[ f_{4321}=d_{43x}\hspace{2mm}d_{21x}+d_{43y}\hspace{2mm}d_{21y}+d_{43z}\hspace{2mm}d_{21z} \]
\[ f_{1321}=d_{13x}\hspace{2mm}d_{21x}+d_{13y}\hspace{2mm}d_{21y}+d_{13z}\hspace{2mm}d_{21z} \]
\[ f_{4343}=d_{43x}\hspace{2mm}d_{43x}+d_{43y}\hspace{2mm}d_{43y}+d_{43z}\hspace{2mm}d_{43z} \]
\[ f_{2121}=d_{21x}\hspace{2mm}d_{21x}+d_{21y}\hspace{2mm}d_{21y}+d_{21z}\hspace{2mm}d_{21z} \]

Con questi cinque fattori si calcolano altri due valori che costituiranno il numeratore e il denominatore di una successiva frazione:

\[ numer = f_{1343}\hspace{2mm}f_{4321}-f_{1321}\hspace{2mm}f_{4343} \]
\[ denom = f_{2121}\hspace{2mm}f_{4343}-f_{4321}\hspace{2mm}f_{4321} \]

Si calcola poi il coefficiente angolare μa della retta a tra i punti P1P2 quale quoziente della frazione:

\[ μ_a=\frac{numer}{denom} \]

con il quale si ricava anche il coefficiente angolare μb della retta b tra i punti P3P4:

\[ μ_b=\frac{f_{1343}+f_{4321}\hspace{2mm}μ_a}{f_{4343}} \]

Con i coefficienti angolari delle due rette si trovano le coordinate dei due estremi Pa e Pb del minimo segmento (perpendicolare a entrambe) che le connette:

\[ x_a=x_1+μ_a\hspace{2mm}d_{21x}\hspace{10mm}y_a=y_1+μ_a\hspace{2mm}d_{21y}\hspace{10mm}z_a=z_1+μ_a\hspace{2mm}d_{21z} \]
\[ x_b=x_3+μ_b\hspace{2mm}d_{43x}\hspace{10mm}y_b=y_3+μ_b\hspace{2mm}d_{43y}\hspace{10mm}z_b=z_3+μ_b\hspace{2mm}d_{43z} \]

Con le coordinate di questi due punti è opportuno calcolare la lunghezza Dab della proiezione sul piano orizzontale del segmento 3D PaPb così da verificare eventuali errori di rilevazione:

\[ D_{ab}=\sqrt{(x_b-x_a )^2+(y_b-y_a )^2} \]

Se questa lunghezza rientra nella tolleranza imposta, possiamo considerare concluso il calcolo, nel senso che i punti Pa e Pb costituiscono l’iperdeterminazione del punto di intersezione cercato al quale il successivo calcolo del rilievo attribuirà le coordinate definitive con i criteri della compensazione dei rilievi, criteri per i quali rimando i più appassionati alla topografia alla visione del video Introduzione alla Teoria degli Errori in Topografia di TopGeometri, al punto in cui viene trattata la compensazione di un rilievo topografico.

Bene, descritto l’algoritmo, può essere come al solito utile vederlo applicato numericamente con un esempio concreto. Lo farò qui di seguito con riferimento al rilievo riprodotto in Figura 6.

Figura 6L’intersezione 3D nelle tabelle di Geocat: le stazioni 100 e 200 sono state fatte su punti GPS e orientate su altrettanti punti GPS (1001 e 1002), dopodiché da ciascuna è stato osservato solo angolarmente il punto 101

Il primo passaggio è il calcolo della distanza tra le stazioni 100 e 200 applicando il teorema di Pitagora alle coordinate determinate dal GPS (sono desumibili dalla tabella coordinate in basso di Figura 6):

\[ D_{100-200}=\sqrt{(x_{200}-x_{100})^2+(y_{200}-y_{100})^2 } \]
\[ D_{100-200}=\sqrt{(-1653.154-312.857)^2+(215.457-296.657)^2}=1967.687 \]

Poi dobbiamo calcolare la quota  del centro strumentale di ciascuna delle due stazioni sommando alla quota del punto a terra (determinato dal calcolo GPS) la rispettiva altezza strumentale (detto anche “piano di mira”, evidenziata nei due riquadri in rosso a sinistra della tabella celerimetrica di Figura 6):

\[ z_{100\hspace{1mm}st}= z_{100\hspace{1mm}terra}+hs_{100}\hspace{20mm}z_{200\hspace{1mm}st}= z_{200\hspace{1mm}terra}+hs_{200} \]
\[ z_{100\hspace{1mm}st}=27.269+1.480=28.749\hspace{10mm}z_{200\hspace{1mm}st}=193.959+1.580=195.539 \]

Ora calcoliamo le coordinate planimetriche x e y dei due punti fittizi P2 e P4 che ci servono per definire le due rette 3D. Troviamo dapprima l’azimut tra ciascuna stazione e il punto di intersezione incognito, sottraendo dall’angolo orizzontale letto in campagna la relativa correzione angolare calcolata dall’elaborazione dei rilievi GPS e TS ed evidenziata nei due riquadri in rosso a destra della tabella celerimetrica di Figura 6:

\[ ϑ_{100-101}=α_{100-101}-corr_{100}\hspace{20mm}ϑ_{200-101}=α_{200-101}-corr_{200} \]
\[ ϑ_{100-101}=82.8273-106.6802=-23.8529+400=376.1470 \]
\[ ϑ_{200-101}=172.4941-97.2535=75.2402 \]

Possiamo ora calcolare le coordinate x, y, z dei punti fittizi P2 e P4:

\[ x_2=x_{100}+D_{100-200}\hspace{3mm}sinϑ_{100-101} \]
\[ y_2=y_{100}+D_{100-200}\hspace{3mm}cosϑ_{100-101} \]
\[ z_2=x_{100\hspace{1mm}st}+D_{100-200}\hspace{3mm}cosφ_{100-101} \]
\[ x_2=312.857+1967.687\hspace{3mm}sin⁡\hspace{1mm}376.1470=-407.271 \]
\[ y_2=296.657+1967.687\hspace{3mm}cos\hspace{1mm}⁡376.1470=2127.834 \]
\[ z_2=28.749+1967.687\hspace{3mm}cos\hspace{1mm}⁡99.3471=48.929 \]
\[ x_4=x_{200}+D_{100-200}\hspace{3mm}sinϑ_{200-101} \]
\[ y_4=y_{200}+D_{100-200}\hspace{3mm}cosϑ_{200-101} \]
\[ z_4=x_{200\hspace{1mm}st}+D_{100-200}\hspace{3mm}cosφ_{200-101} \]
\[ x_4=-1653.154+1967.687\hspace{3mm}sin⁡\hspace{1mm}75.2402=167.580 \]
\[ y_4=215.457+1967.687\hspace{3mm}cos⁡\hspace{1mm}75.2402=961.593 \]
\[ z_4=195.539+1967.687\hspace{3mm}cos\hspace{1mm}⁡105.4784=26.420 \]

Troviamo ora le differenze di coordinate tra i punti:

\[ d_{13}=P_3-P_1\hspace{10mm}d_{21}=P_1-P_2\hspace{10mm}d_{43}=P_3-P_4 \]

dove P1 e P3 sono le due stazioni 100 e 200, mentre P2 e P4 sono i due punti fittizi appena calcolati:

\[ d_{13x}=x_{100}-x_{200}=312.857-(-1653.154)=1966.011 \]
\[ d_{13y}=y_{100}-y_{200}=296.657-215.457=81.200 \]
\[ d_{13z}=z_{100}-z_{200}=28.749-195.539=-166.790 \]
\[ d_{21x}=x_2-x_{100}=-407.271-312.857=-720.127 \]
\[ d_{21y}=y_2-y_{100}=2127.834-296.657=1831.177 \]
\[ d_{21z}=z_2-z_{100}=48.929-28.749=20.180 \]
\[ d_{43x}=x_4-x_{200}=167.580-(-1653.154)=1820.734 \]
\[ d_{43y}=y_4-y_{200}=961.593-215.457=746.136 \]
\[ d_{43z}=z_4-z_{200}=26.420-195.539=-169.119 \]

Con questi valori calcoliamo i cinque fattori f:

\[ f_{1343} = d_{13x}\hspace{3mm}d_{43x}+d_{13y}\hspace{3mm}d_{43y}+d_{13z}\hspace{3mm}d_{43z} \]
\[ f_{1343} = 1966.011∙1820.734 + 81.200∙ 746.136 + -166.790 ∙ -169.119 \]
\[ f_{1343}=3668376.567 \]
\[ f_{4321}=d_{43x}\hspace{3mm}d_{21x}+d_{43y}\hspace{3mm}d_{21y}+d_{43z}\hspace{3mm}d_{21z} \]
\[ f_{4321}=1820.734 ∙ -720.127+746.136 ∙ 1831.177+ -169.119 ∙ 20.180 \]
\[ f_{4321}=51733.094 \]
\[ f_{1321}=d_{13x}\hspace{3mm}d_{21x}+d_{13y}\hspace{3mm}d_{21y}+d_{13z}\hspace{3mm}d_{21z} \]
\[ f_{1321}=1966.011∙-720.127 + 81.200∙ 1831.177 + -166.790 ∙ 20.180 \]
\[ f_{1321}=-1270451.540 \]
\[ f_{4343}=d_{43x}\hspace{3mm}d_{43x}+d_{43y}\hspace{3mm}d_{43y}+d_{43z}\hspace{2mm}d_{43z} \]
\[ f_{4343}=1820.734∙1820.734 + 746.136∙ 746.136 +-169.119 ∙ -169.119 \]
\[ f_{4343}=3900392.629 \]
\[ f_{2121}=d_{21x}\hspace{3mm}d_{21x}+d_{21y}\hspace{3mm}d_{21y}+d_{21z}\hspace{2mm}d_{21z} \]
\[ f_{2121}=-720.127∙-720.127 + 1831.177∙ 1831.177 + 20.180 ∙ 20.180 \]
\[ f_{2121}=3872198.478 \]

E con questi cinque fattori calcoliamo il numeratore e il denominatore della successiva frazione:

\[ numer=f_{1343}\hspace{3mm}f_{4321}-f_{1321}\hspace{3mm}f_{4343} \]
\[ numer=3668376.567 ∙ 51733.094\hspace{3mm}-\hspace{3mm}-1270451.540 ∙ 3900392.629 \]
\[ numer=5145036289946.85 \]
\[ denom=f_{2121}\hspace{3mm}f_{4343}-f_{4321}\hspace{3mm}f_{4321} \]
\[ denom=3872198.478 ∙ 3900392.629\hspace{3mm}-\hspace{3mm}51733.094 ∙ 51733.094 \]
\[ denom=15100418088556.90 \]

Quindi calcoliamo il coefficiente angolare μa della retta a (P1P2), cioè il raggio di mira uscente dalla stazione 100, quale quoziente della frazione:

\[ μ_a=\frac{5145036289946.85}{15100418088556.90}=0.340721446239013 \]

E poi il coefficiente angolare μb della retta b (P3P4), cioè il raggio di mira uscente dalla stazione 200:

\[ μ_b=\frac{f_{1343}+f_{4321}\hspace{3mm}μ_a}{f_{4343}} \]
\[ μ_b=\frac{3668376.567+51733.094∙0.340721446239013}{3900392.629} \]
\[ μ_b=0.945033870183323 \]

Con i coefficienti angolari delle due rette troviamo finalmente le coordinate dei due estremi Pa e Pb del minimo segmento che le connette:

\[ x_a=x_{100}+μ_a\hspace{3mm}d_{21x} \]
\[ x_a=312.857+0.340721446239013 ∙ -720.127=67.494 \]
\[ y_a=y_{100}+μ_a\hspace{3mm}d_{21y} \]
\[ y_a=296.657+0.340721446239013∙1831.177=920.578 \]
\[ z_a=z_{100}+μ_a\hspace{3mm}d_{21z} \]
\[ z_a=28.749+0.340721446239013 ∙ 20.180=35.624 \]
\[ x_b=x_{200}+μ_b\hspace{3mm}d_{43x} \]
\[ x_b=-1653.154+0.945033870183323 ∙ 1820.734=67.502 \]
\[ y_b=y_{200}+μ_b\hspace{3mm}d_{43y} \]
\[ y_b=215.457+0.945033870183323 ∙ 746.136=920.580 \]
\[ z_b=z_{200}+μ_b\hspace{3mm}d_{43z} \]
\[ z_b=195.539+0.945033870183323 ∙ -169.119=35.715 \]

A verifica sulla bontà delle rilevazioni calcoliamo la distanza piana Dab del segmento Pa-Pb:

\[ D_{ab}=\sqrt{(x_b-x_a)^2+(y_b-y_a)^2} \]
\[ D_{ab}=\sqrt{(67.502-67.494)^2+(920.580-920.578)^2}=0.008 \]

Questo risultato ci conferma la correttezza sia delle osservazioni eseguite in campagna che quella dei calcoli sopra esposti e ci consente di mantenere le coordinate di entrambi i punti Pa e Pb per il calcolo complessivo finale di compensazione dell’intero rilievo.

[1] Che poi, se devo essere sincero, io trovo molto di più nei miei vecchi testi e quaderni di topofrafia che ho gelosamente conservato dalla scuola di geometra.

%d blogger hanno fatto clic su Mi Piace per questo: