[php] trilaterazione, localizzazione sorgente GSM, WIFI, RADIO, etc..

[Rubberick]

ciao a tutti :)

mi sto divertendo a giochicciare con google earth e i suoi kml...

ho a disposizione una serie di dati sulle celle 2g/3g della zona, ho gia' scriptato tutto quello che mi serve tranne una funzione di trilaterazione per la localizzazione + accurata delle torrette

in pratica ho una lista di valori

latitudine, longitudine, altezza (in m), potenza (in dBm), accuratezza gps (in m)

dovrei realizzare una funzione che utilizzando un minimo di 3 o 4 record a salire mi restituisca principalmente latitudine e longitudine (e possibilmente in un secondo momento altezza) della sorgente di segnale basandomi sulle misurazioni

andrei per step...

prima di tutto realizzare una trilaterazione solo latitudine / longitudine basandomi su 3 punti

poi lo estendo su + punti realizzando una MULTILATERAZIONE

poi faccio un discorso di pesi, dove do maggior peso ai record che mi danno maggiore potenza in dBm e discorso inverso per i punti che mi danno minore accuratezza (ad es 50m invece che 3m)

ho trovato uno snippet di codice in phyton e stavo cercando di convertirlo in php ma non ho ben capito tutti i passaggi... qualcuno mi aiuta a capirli un po' meglio? inoltre usa un riferimento sferico se ho capito e non ellissoide

non so bene cosa è numpy, ma comunque pare che si debbano fare un pò di calcoli con matrici... in matlab credo il problema sarebbe realizzabile in pochi min... in php diventa un pò + complesso...

from math import *
from numpy import *

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
# 1. Convert Lat/Long to radians
# 2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = array([xA, yA, zA])
P2 = array([xB, yB, zB])
P3 = array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i*x)/j)

# only one case shown here
z = sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon`

[Rubberick]

ora che ci penso ma è corretto questo ragionamento?

che sono DistA DistB e DistC?

io ho delle rilevazioni basate sul mio gps quindi la posizione sicura (a meno dell'accuratezza del gps) e' la mia non quella della torre ^^ altrimenti avrei risolto

io alla fine ho gia' tutti i punti che mi servono...

potrei fare una media di lat e lon pesata su potenza e accuratezza?

viene uno schifo eh?

[Rubberick]

ok, credo di avere capito sono o i raggi o i diametri delle sorgenti...

che a questo punto potrebbero essere basati direttamente sulla potenza tramite la perdita in free loss dal segnale

[Rubberick]

uff mi sto incazzando come una iena complice anche il caldo per fare un calcolo quantomeno banale...

sto cercando di fare una funzioncina semplice x il calcolo della distanza in base ai decibel..

è un primo step per iniziare poi la raffino considerando l'attenuazione data da possibili ostacoli

In spazio libero la faccenda dovrebbe funzionare così:

http://en.wikipedia.org/wiki/Free-space_path_loss

la formula e' la seguente:

http://upload.wikimedia.org/math/a/d/4/ad4e48cb0ee25df10f76e46171df3974.png

con d espressa in metri e f in Hz

benissimo pertanto faccio la formula inversa

dB = 20log10(d) + 20log10(f) - 147.55

dB/20 = log10(d) + log10(f) - 147.55/20

dB/20 = log10(d) + log10(f) - 7.3775

passando i 2 termini dall'altra parte...

dB/20 - log10(f) + 7.3775 = log10(d)

fin qui giusto no? sto facendo un passaggino alla volta x non sbagliare

a questo punto semplicemente faccio 10^ da un lato e dall'altro

10^( dB/20 - log10(f) + 7.3775 ) = 10^( log10(d) )

10^( dB/20 - log10(f) + 7.3775 ) = d

dove ricordiamo la freq e' espressa in Hz e la d in metri

quindi mi preparo una funzioncina

db_distance($db,$f) {
return pow(10,($db/20 + 7.3775 - log10($f)));
}

da qui se faccio un

echo db_distance(-80,2.45*pow(10,9));

(ho usato 2.45 GHz per provare il caso WIFI)

mi esce 0.000000973495426720451714742120 m

ora...

se prendo una funzione già fatta da un tizio che e' valida solo per il wifi (2.45GHz)

function db_km($db){
$cst=-40.4;
return pow(10,($cst-$db)/20)/1000;
}

mi da 0.095499258602144 km

che in m sono 0.000095499258602144 m

mi trovo 2 ordini di scarto... in cosa ho sbagliato? O_o'

prima di usarla x le altre frequenze 850MHz, 900Mhz, 1800Mhz etc.. vorrei essere sicuro che la funzione sia corretta..

[Rubberick]

forse ho trovato il problema c'e' un doppio meno che va considerato per l'attenuazione...

non si puo' inserire -80 db ma tocca usare 80 o altrimenti se si vuole usare i -80 va messo un - davanti a $db


function db_distance($db,$f) {
return pow(10,(-$db/20 + 7.3775 - log10($f)));
}

e cosi' funziona correttamente :)

lo so sto facendo un monologo ma prima o poi qualcuno interessato comparira' xD

[Luciano69]

Ciao
posso sapere le tue info di partenza
le cordinate delle celle GSM o UMTS come le hai trovate?