Вычисление координат точки

Доброе время суток.

Столкнулся с такой задачей: есть координаты объекта с приемника GPS. На нем установлен электронный компас, датчики положения (угол крена и тангажа) и дальномер. Это устройство смотрит на какую-то точку Б. Необходимо по имеющимся данным вычислить координаты точки Б.

К сожалению к геодезии не имею никакого отношения, поэтому не знаю куда копать. Пожалуйста, подскажите по каким формулам можно посчитать?

 

Ну если учесть, что GPS/Глонасс дает приличную ошибку, наш компас еще больше и т.д. То думаю хорошая точность не получится все равно. Поэтому можно считать, что сойдет низкая. Конкретно +- сколько метров сказать не могу. Но чтобы полученная координата имела какое-то отношение к действительности )

 

осмелюсь предположить, что найти участок на местности... Или проверить докУменты...

Это второе. Первое - в какой стистеме координат...

 

Это макет устройства обнаружения и сопровождения целей, обнаружения несанкционированного проникновения на охраняемую территорию и т.д. Координаты цели, хотя бы примерно, нужно указывать. Координаты- широта, долгота, как с навигатора.

 

Нет, только точки.

Есть собственные координаты. Есть направление, куда смотрит прибор. Есть расстояние до объекта, на который он смотрит. Есть углы наклона прибора. Нужна координата той точки, куда он смотрит.

 

Я чуть поправлю Шуфотинского. Приращения к прямоугольным координатам несколько сложно "прилепить" к исходным геодезическим координатам (широта, долгота). Здесь надо составить простую таблицу: сколько секунд (минут) по широте или долготе приходится на 1, 10, 100 и т.д. метров, приращения прямоугольных координат переводить в секунды, и уже их прибавлять к исходным координатам.

 

Ну вот он показывает- 5543.3380N 03733.6810E Какая это система координат?

 

Посмотрите руководство пользователя вашего навигатора и скорее всего это будет в системе координат WGS-84 (если у вас в навигаторе нет возможности пользователю изменять систему координат), что будет означать: 55 градусов 43,338 минуты северной широты и 37 градусов 33,681 восточной долготы.

 

Спасибо! Кажется разобрался. С приближениями, конечно. Но ошибка в расчетах получается приемлемой (проверил по гугл-картам). Теперь осталось заложить алгоритм в девайс...

 

Успехов!

 

Ну как обычно, геодезисты все усложнили, затеяли спор о 2мм. Начали рассказывать о проекциях и формы Земли.
Понимаю, что речь идет о 10 летнем посте. И скорее всего задача у того парня решена.
Но может кому то еще понадобиться. Речь скорее всего шла о каком то локаторе, на выходе надо получить координаты точки в WGS84/
Скорее всего наблюдаемый аппарат фиксирован, что упрощает задачу.
Это все сводиться к прямой геодезической задачи. в которой мы будем искать ΔШироты и ΔДолготы
Для начала определяем соотношение сколько в одном метре у нас градусов по широте и по долготе в конкретной точки (области где производим работы)
Их можно принять за константы.
В начале все решаем в плоскости... находим Δx Δy далее переводим в угловые приращения. Далее прибавляем исходным координатам.
Узнать соотношение км/градусы в долготе и в широте конкретно в нужной точки (области) можно из множество калькуляторов в интернете. Задав поочередно к примеру. 45 00 00 с ш. 50 00 00 в.д (исходная точка) ищем расстояние по долготе вводим точку 46 00 00 с ш. 50 00 00 в.д получили расстояние поделил, и далее по долготе. 45 00 00 с ш. 51 00 00

Если вычислительные мощности позволяют, можно переводить в проекцию и обратно...
Либо загорячиться, на формулу Винсента, она гуглица за раз.
Либо воспользоваться моей. На яве. /Весь код рабочий, точность высокая. Сейчас уже не скажу... надо смотреть но вопрос первых 0.3мм, и сходимость до 3 мм если крутить прямая обратная прямая.

//расстояние на элепсоиде через формулу Венсенти
// где B широта L долгота в градуса, десятичных двух точек
// Aa большая полуось э-да Ea, Ea сжатие
// на выходе расстояние.
// Обратная геодезическая задача
public static double vincenty(double B, double L, double B1,
double L1, double Aa, double Ea) {
// Параметры эллипсоида:
double a = Aa;//6378245.0;
double f = Ea;//1 / 298.3;
double b = (1 - f) * a;
double EPS = 0.5*Math.pow(10,-30);//0.5E-30;
double APARAM;
double BPARAM, CPARAM, OMEGA, TanU1, TanU2, lambda, LambdaPrev, SinL, CosL, USQR, U1, U2, SinU1, CosU1, SinU2, CosU2, SinSQSigma, CosSigma, TanSigma, Sigma, SinAlpha, Cos2SigmaM, DSigma;
double PI = Math.PI;
L = L * (PI / 180);
B = B * (PI / 180);
L1 = L1 * (PI / 180);
B1 = B1 * (PI / 180); // Пересчет значений координат в
// радианы
TanU1 = (1 - f) * Math.tan(B);
TanU2 = (1 - f) * Math.tan(B1);
U1 = Math.atan(TanU1);
U2 = Math.atan(TanU2);
SinU1 = Math.sin(U1);
CosU1 = Math.cos(U1);
SinU2 = Math.sin(U2);
CosU2 = Math.cos(U2);
OMEGA = L1 - L;
lambda = OMEGA;
do { // Начало цикла итерации
LambdaPrev = lambda;
SinL = Math.sin(lambda);
CosL = Math.cos(lambda);
SinSQSigma = (CosU2 * SinL * CosU2 * SinL)
+ (CosU1 * SinU2 - SinU1 * CosU2 * CosL)
* (CosU1 * SinU2 - SinU1 * CosU2 * CosL);
CosSigma = SinU1 * SinU2 + CosU1 * CosU2 * CosL;
TanSigma = Math.sqrt(SinSQSigma) / CosSigma;
if (TanSigma > 0) {
Sigma = Math.atan(TanSigma);
} else {
Sigma = Math.atan(TanSigma) + PI;
}
if (SinSQSigma == 0) {
SinAlpha = 0;
} else {
SinAlpha = CosU1 * CosU2 * SinL / Math.sqrt(SinSQSigma);
}
if ((Math.cos(Math.asin(SinAlpha)) * Math.cos(Math.asin(SinAlpha)) == 0)) {
Cos2SigmaM = 0;
} else {
Cos2SigmaM = CosSigma
- (2 * SinU1 * SinU2 / (Math.cos(Math.asin(SinAlpha)) * Math
.cos(Math.asin(SinAlpha))));
}
CPARAM = (f / 16.0)
* Math.cos(Math.asin(SinAlpha))
* Math.cos(Math.asin(SinAlpha))
* (4 + f
* (4 - 3 * Math.cos(Math.asin(SinAlpha))
* Math.cos(Math.asin(SinAlpha))));
lambda = OMEGA
+ (1 - CPARAM)
* f
* SinAlpha
* (Math.acos(CosSigma) + CPARAM
* Math.sin(Math.acos(CosSigma))
* (Cos2SigmaM + CPARAM * CosSigma
* (-1 + 2 * Cos2SigmaM * Cos2SigmaM)));
} while ((Math.sqrt((lambda - LambdaPrev) * (lambda - LambdaPrev))) > EPS);
// Конец цикла итерации
USQR = Math.cos(Math.asin(SinAlpha)) * Math.cos(Math.asin(SinAlpha))
* (a * a - b * b) / (b * b);
APARAM = 1 + (USQR / 16384)
* (4096 + USQR * (-768 + USQR * (320 - 175 * USQR)));
BPARAM = (USQR / 1024)
* (256 + USQR * (-128 + USQR * (74 - 47 * USQR)));
DSigma = BPARAM
* Math.sqrt(SinSQSigma)
* (Cos2SigmaM + BPARAM
/ 4
* (CosSigma * (-1 + 2 * Cos2SigmaM * Cos2SigmaM) - BPARAM
/ 6
* Cos2SigmaM
* (-3 + 4 * SinSQSigma)
* (-3 + 4 * Cos2SigmaM * Cos2SigmaM)));
return b * APARAM * (Sigma - DSigma);
}
// Венсете дубль два Обратная задача. http://www.planetcalc.ru/722/
// На входе долгота широта в градусах, Дистанция, угол выхода, угол входа
// Обратная геодезическая задча
public static double[] Vincenty (double la1,double lo1,double la2,double lo2, double Aa,double Ea ){

// так.
// https://en.wikipedia.org/wiki/Vincenty's_formulae
// Расчет между двумя точками на земной поверхности
double a=Aa; //length of semi-major axis of the ellipsoid (radius at equator); (6378137.0 metres in WGS-84)
double f=Ea; //ƒ flattening of the ellipsoid; (1/298.257223563 in WGS-84)
//double b = (1 - f)*a; //length of semi-minor axis of the ellipsoid (radius at the poles); (6356752.314245 meters in WGS-84)
double Φ1=Math.toRadians(la1); // latitude of the points 1;
double Φ2=Math.toRadians(la2); // latitude of the points 1;
double U1 = Math.atan((1-f)*Math.tan (Φ1));// reduced latitude (latitude on the auxiliary sphere)
double U2 = Math.atan((1-f)*Math.tan (Φ2));// reduced latitude (latitude on the auxiliary sphere)
double L1=Math.toRadians (lo1);
double L2=Math.toRadians (lo2);
double L = L2 - L1; // difference in longitude of two points (разница в долготе двух точек)
double s;// ellipsoidal distance between the two points;
// f=(geoid.max-geoid.min)/geoid.max;
// f*geoid.max=geoid.max-geoid.min
// geoid.max-f*geoid.max=geoid.min

double f_16=f/16;
double Bb=Aa-Aa*f;// Малая полуось
double b2 = Bb*Bb;
double e_inv = (a*a-b2)/b2;
double minDelta = 0.000000000000000001;
U1 = Math.atan((1-f)*Math.tan(Math.toRadians(la1)));
U2 = Math.atan((1-f)*Math.tan(Math.toRadians(la2)));
double sinU1 = Math.sin(U1);
double sinU2 = Math.sin(U2);
double cosU1 = Math.cos(U1);
double cosU2 = Math.cos(U2);
L = Math.toRadians(lo2-lo1);
double cos2SigmaM2;
double lambda = L;
double dLambda = 0;
double sinAlpha;
double cos2SigmaM;
while(true) {
double sinLambda = Math.sin(lambda);
double cosLambda = Math.cos(lambda);
double sin2Sigma = Math.pow(cosU2*sinLambda,2)+Math.pow(cosU1*sinU2-sinU1*cosU2*cosLambda,2);
double sinSigma = Math.sqrt(sin2Sigma);
double cosSigma = sinU1*sinU2+cosU1*cosU2*cosLambda;
double sigma = Math.atan2(sinSigma,cosSigma);

if (sinSigma !=0){ sinAlpha=cosU1*cosU2*sinLambda/sinSigma;}//проверка дел на ноль
else{sinAlpha=0;};

double cos2a = 1-sinAlpha*sinAlpha;

if (cos2a !=0){ cos2SigmaM =cosSigma - 2*sinU1*sinU2/cos2a;}//проверка дел на ноль
else{cos2SigmaM=0;};

double C = f_16*cos2a*(4+f*(4-3*cos2a));
double oldLambda = lambda;
cos2SigmaM2 = cos2SigmaM*cos2SigmaM;
lambda = L+(1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM2)));
dLambda = lambda-oldLambda;
if ( Math.abs(dLambda)<=minDelta ) {
double u2 = cos2a*e_inv;
double A =1+ (u2/16384)*(4096+u2*(-768+u2*(320-175*u2)));
double B =(u2/1024)*(256+u2*(-128+u2*(74-47*u2)));
double dSigma = B*sinSigma*(cos2SigmaM+(B/4)*(cosSigma*(-1+2*cos2SigmaM2)-(B/6)*cos2SigmaM*(-3+4*sin2Sigma)*(-3+4*cos2SigmaM2)));
s = Bb*A*(sigma-dSigma);// судя запихнють правильный малый радиус земли
cosLambda = Math.cos(lambda);
sinLambda = Math.sin(lambda);
double angle1 = Math.toDegrees(Math.atan2( cosU2*sinLambda,cosU1*sinU2-sinU1*cosU2*cosLambda));
double angle2 = Math.toDegrees(Math.atan2( cosU1*sinLambda,-sinU1*cosU2 + cosU1*sinU2*cosLambda));

double []r = {s,angle1,angle2};
return r;}
}
}
// По координате угле и расстоянию находим вторую точку на геоиде в географических координатах.
// прямая геодезическая задача на элепсоиде
public static double[] VincentyII (double la1, double lo1, double a1,
double s, double Aa, double Ea) {
// var ellipsoids = { "wgs84" : { "min":6356752.314, "max":6378137 }, "sk42" : { "min":6356863, "max":6378245 }, "sphere" : { "min":6378137, "max":6378137 } };
// var geoid = ellipsoids[ellipsoid];
double f = Ea; // (geoid.max-geoid.min)/geoid.max;
double f_16=f/16;
// малая полуось b = (1 - ƒ) a
double b = (1 - f) * Aa;
double b2 = b*b;
double e_inv = (Aa*Aa-b2)/b2;
// Параметры эллипсоида:
//double a = Aa;//6378245.0;
f = Ea;//1 / 298.3;
// double b = (1 - f) * a;
double EPS = 0.5*Math.pow(10,-30);//0.5E-30;
a1=Math.toRadians(a1);// точно ли угол в Радианах?
double U1 = Math.atan((1-f)*Math.tan(Math.toRadians((la1))));
double sinU1 = Math.sin(U1);
double cosU1 = Math.cos(U1);
double sinAlpha = cosU1*Math.sin(a1);
double sinAlpha2 = sinAlpha*sinAlpha;
double cosAlpha2 = 1-sinAlpha2;
double u_2 = cosAlpha2*e_inv;
double A = 1+u_2/16384*(4096+u_2*(-768+u_2*(320-175*u_2)));
double B = u_2/1024*(256+u_2*(-128+u_2*(74-47*u_2)));
double sigma = s/(b*A);
double sigmaBase = sigma;
double cosa1 = Math.cos(a1);
double sina1 = Math.sin(a1);
double sigma1 = Math.atan2(Math.tan(U1),cosa1);

while(true) {
double twoSigmaM = 2*sigma1+sigma;
double cos2SigmaM = Math.cos(twoSigmaM);
double cos2SigmaM2 = cos2SigmaM*cos2SigmaM;
double sinSigma = Math.sin(sigma);
double cosSigma = Math.cos(sigma);
double deltaSigma = B*sinSigma*( cos2SigmaM+1/4*B*(cosSigma*(-1+2*cos2SigmaM2)-1/6*B*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM2) ) );
double oldSigma = sigma;
sigma = sigmaBase+deltaSigma;
double dSigma = Math.abs(sigma-oldSigma);

if ( dSigma<EPS) {
double [] res = {0,0,0};
res[0] = Math.toDegrees(Math.atan2( sinU1*cosSigma+cosU1*sinSigma*cosa1,(1-f)*Math.sqrt(sinAlpha2+Math.pow((sinU1*sinSigma-cosU1*cosSigma*cosa1),2)) ));
double lambda = Math.atan2( sinSigma*sina1,cosU1*cosSigma-sinU1*sinSigma*cosa1);
double C = f_16 * cosAlpha2 * (4 + f * (4 - 3 * cosAlpha2));
double L = lambda-(1-C)*f*sinAlpha*(sigma +C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+ 2*cos2SigmaM2)));
res[1] = lo1 + Math.toDegrees(L);
res[2] = Math.atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosa1);
return res;
}
}
}





// прямая геодезическая задача на элепсоиде
// По координате угле и расстоянию находим вторую точку на геоиде в географических координатах.
public static double[] vincentyII ( double la1, double lo1, double a1, double s, double Aa, double Ea) {
double minDelta = 0.000000000000000001;
a1=Math.toRadians(a1);
double f_16=Ea/16;
double Bb=Aa-Aa*Ea;// Малая полуось
double b2 = Bb*Bb;
double U1 = Math.atan((1-Ea)*Math.tan(Math.toRadians((la1))));
double sinU1 = Math.sin(U1);
double cosU1 = Math.cos(U1);
double sinAlpha = cosU1*Math.sin(a1);
double sinAlpha2 = sinAlpha*sinAlpha;
double cosAlpha2 = 1-sinAlpha2;
double e_inv = (Aa*Aa-b2)/b2;

double u_2 = cosAlpha2*e_inv;
double A = 1+u_2/16384*(4096+u_2*(-768+u_2*(320-175*u_2)));
double B = u_2/1024*(256+u_2*(-128+u_2*(74-47*u_2)));
double sigma = s/(Bb*A);
double sigmaBase = sigma;
double cosa1 = Math.cos(a1);
double sina1 = Math.sin(a1);
double sigma1 = Math.atan2(Math.tan(U1),cosa1);
while(true) {
double twoSigmaM = 2*sigma1+sigma;
double cos2SigmaM = Math.cos(twoSigmaM);
double cos2SigmaM2 = cos2SigmaM*cos2SigmaM;
double sinSigma = Math.sin(sigma);
double cosSigma = Math.cos(sigma);
double deltaSigma = B*sinSigma*( cos2SigmaM+1/4*B*(cosSigma*(-1+2*cos2SigmaM2)-1/6*B*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM2) ) );
double oldSigma = sigma;
sigma = sigmaBase+deltaSigma;
double dSigma = Math.abs(sigma-oldSigma);
if ( dSigma<minDelta) {
double []res={0,0,0} ;
res[0] = Math.toDegrees(Math.atan2( sinU1*cosSigma+cosU1*sinSigma*cosa1,(1-Ea)*Math.sqrt(sinAlpha2+Math.pow((sinU1*sinSigma-cosU1*cosSigma*cosa1),2)) ));
double lambda = Math.atan2( sinSigma*sina1,cosU1*cosSigma-sinU1*sinSigma*cosa1);
double C = f_16 * cosAlpha2 * (4 + Ea * (4 - 3 * cosAlpha2));
double L = lambda-(1-C)*Ea*sinAlpha*(sigma +C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+ 2*cos2SigmaM2)));
res[1] = lo1 + Math.toDegrees(L);
res[2] =Math.toDegrees( Math.atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosa1));
return res;
}
}
}

 

Stepan_S, ПЕРЕВОД КООРДИНАТ ИЗ СТАРТОВОЙ СИСТЕМЫ В ГЕОДЕЗИЧЕСКУЮ

Хотя и разработаны формулы Висенти были для уменьшения длины кода программы, но использовать их сейчас, при наличии GeographicLib, которая портирована на многие современные языки программирования и является частью PROJ не самая лучшая идея.
Кроме того, есть достаточно элегантное решение (прямой и обратной геодезических задач) 2005 года.

 

вот и приходят ко мне студенты, которые в душе не понимают, что происходит в чуда коробках (тахеометрах, GPSсках) вот и рассказывают, мне "мы же координаты ЖПСКОЙ определили" А когда я говорю, что база - ровер, не определяет координаты, а определяет только вектор... смотрят на меня как бараны на новые ворота. И на вопрос, почему такой результат? "Ну мы библиотеку применили, а что в ней происходит в душе не понимаем" И приносят статику, на которую без слез не взглянешь...
И начинаешь, по сотому кругу рассказывать, что и как и что косинусы маленькие, и про тангенсы 90 градусов...
Элегантное решение надо по изучать, что там нового. Не видел его.

 

вдогонку про 90 градусов из жизни анекдот... ну рассказываю я, что тангенс (это людям после высшей математике, геодезисты) 90 градусов, это не очень хорошая идея... и рассказываю почем, так получилось, типа делить на ноль нельзя по этому и тангенс.... и тут вопрос от "гения" " А почему нельзя делить, то на ноль" .. но я дотошный... решил спросить, "ну дели" (ожидал ответ ∝), но услышал ответ НОЛЬ... //ЗАНАВЕС.