расчет координат точки на местности используя исходные данные модуля навигатора

GAMP: An open-source software of multi-GNSS precise point positioning using undifferenced and uncombined observations

PPPH: A MATLAB-based software for multi-GNSS precise point positioning analysis

goGPS: open-source MATLAB software -- by Antonio M. Herrera, Hendy F. Suhandri, Eugenio Realini, Mirko Reguzzoni, and M. Clara de Lacy

 

Доброе утро. Что то не сходится у меня расчет в это коде. Запускаю с теми исходными данными что заданы. результат сходится. считываю навигационное сообщение на этот же спутник и на это же время но на следующий день. результат сходится. Ввожу данные и время другого спутника, - результат не совпадает.
вот 15 спутник - эфемериды:
* 2019 2 12 8 0 0.00000000
PG15 323.842837 -24393.481094 -10268.435720 -328.217094 8 9 5 84

навигационное сообщение
15 19 2 12 8 0 0.0-0.328210648149D-03 0.170530256582D-11 0.000000000000D+00
0.160000000000D+02-0.331250000000D+01 0.523093213189D-08-0.124802397881D+01
-0.182539224625D-06 0.112132804934D-01 0.928156077862D-05 0.515363610840D+04
0.201600000000D+06 0.296160578728D-06 0.897675791800D+00 0.726431608200D-07
0.928038015777D+00 0.181906250000D+03 0.763414028836D+00-0.831177437988D-08
-0.243938730327D-09 0.100000000000D+01 0.204000000000D+04 0.000000000000D+00
0.200000000000D+01 0.000000000000D+00-0.107102096081D-07 0.160000000000D+02
0.197370000000D+06 0.000000000000D+00 0.000000000000D+00 0.000000000000D+00

расшифровка:
double Crs = -0.331250000000E+01;
double dn = 0.523093213189E-08;
double M0 = -0.124802397881E+01;
double Cuc = -0.182539224625E-06;
double ec = 0.112132804934E-01;
double Cus = 0.928156077862E-05;
double a = 0.515363610840E+04 * 0.515363610840E+04;
double Toe = 0.201600000000E+06;
double Cic = 0.296160578728E-06;
double W0 = 0.897675791800E+00;
double Cis = 0.726431608200E-07;
double i0 = 0.928038015777E+00;
double Crc = 0.181906250000E+03;
double w = 0.763414028836E+00;
double Wdot = -0.831177437988E-08;
double idot = -0.243938730327E-09;
время T - double T = 28800;
результат
X = 23761731.699 Y = 4267023.204 Z = -11757519.146

 

?
T ≡ 0

 

T = 8 *60 * 60 = 28800


Да, при T - 0 все сразу получилось.

 

Сравниваем с данными в SP3 на момент
* 2019 2 12 8 0 0.00000000
…
PG15 323.842837 -24393.481094 -10268.435720

Берём из навигационного файла последовательно на моменты (выделено жирным)
15 19 2 12 5 59 44.0-0.328223221004D-03 0.170530256582D-11 0.000000000000D+00
…
15 19 2 12 8 0 0.0-0.328210648149D-03 0.170530256582D-11 0.000000000000D+00
…
15 19 2 12 10 0 0.0-0.328198540956D-03 0.170530256582D-11 0.000000000000D+00
Получим:
T = 7216.0000000
X = 323855.110 Y = -24393486.501 Z = -10268423.228
T = 0.0000000
X = 323854.697 Y = -24393486.153 Z = -10268423.669
T = -7200.0000000
X = 323854.624 Y = -24393486.029 Z = -10268423.571

 

Добрый день. Вот время Toe, которое в эфемеридах, - оно там представлено на начало или конец интервала действия навигационного сообщения? вот я смотрю toe на 19 число с 2 до 4 UTC, там 1800000 у меня в это время в 9-45 (новосибирск, разница 7 часов), 187200. (это данные от модуля). Отчего вопрос? У меня время (Врем¤ недели измерени¤ по локальному времени приемника. rcvTow \ в мс ) меньше чем toe. цифры такие, - если toe180000, то Tow 179135, в другом случае, - 194400 - 190784. Мне кажется, так быть не должно.

 

Все таки не понимаю....... везде пишут: t_oe - опорная эпоха (время навигационного сообщения) в секундах от начала недели смотрю часы: вот время на этот момент - 265444s. беру данные от навигатора (в сообщении псевдодальности) примерно в это же время - сообщение Tow 265ххх секунд, Toe из навигатора с сообщении эфемерид - 266400. Но ведь не должна же быть опорная эпоха больше, чем настоящее время?

 

Почему?

--- Сообщения объединены, 20 фев 2019, Оригинальное время сообщения: 20 фев 2019 ---

Воспользуйтесь советом Александр Яковченко,

Обратите внимание на раздел 20.3.4.5 Reference Times. и фразу "Note that a user applying current navigation data will normally be working with negative values of (t-toc) and (t-toe) in evaluating the expansions."

 

1. Эпоха в астрономии (от греч. έποχή — «остановка») — момент времени, для которого определены астрономические координаты или элементы орбиты.
2. t_oe - опорная эпоха (время навигационного сообщения) в секундах от начала недели То есть в моем понимании, - от начала, - до начала следующей эпохи.
другого объяснения не знаю.

--- Сообщения объединены, 20 фев 2019, Оригинальное время сообщения: 20 фев 2019 ---

спасибо. перечитаю. до этого я не дочитал...... значит все правильно пока.........

 

Отличие астрономического термина эпоха от "исторического" (и бытового, без всякого отрицательного смысла) в том, что в астрономии эпоха – это просто конкретный момент времени (любой), а в истории эпоха – интервал времени между моментами. И как момент времени его совсем не обязательно связывать с координатами или любым фазовым вектором состояния любого тела. Поэтому совсем нет смысла говорить о "начале следующей эпохи". Это начало следующей эпохи совсем не нужно, да и его может просто не быть.
P.S. Русская вики по своему обыкновению врёт, чаще в деталях.

 

Добрый вечер. Подскажите пожалуйста, - какой точности вычисления координат спутников достаточно для расчета координат точки на местности, с точностью в пределах плюс/минус 100 метров.

 

https://www.google.ru/search?newwin.......0....1..gws-wiz.......0i22i30.Qpyg7BbONH0

 

Доброе утро. применил поправки.

Код:

  void calk_SV1(int x) // расчет положения спутника по его эфемеридам
  {
    double E;
    a_bp = pow(GP[x].square_big_axis,2); // вычмсляем большую полуось
    n_o = sqrt(gp_earth/(pow(a_bp,3))); //
    n = n_o + GP[x].delta_n; // находим среднее движение
 
   // t = t - ( PR_sv[x] / c );//
    t_k = t - GP[x].t_oe ;   // вычисляем промежуток времени от опорной эпохи эфемерид
     if ( t_k > 302400)       // если t больше 302400 секунд, вычтите 604800 секунд из T.
     t_k = t_k - 604800;
     else if ( t_k < -302400) // Если T меньше -302400 секунд, добавьте 604800 секунд к T.
     t_k = t_k + 604800;
 
    M_sr = GP[x].Mo + (n * t_k); // определяем среднюю аномалию
    E1 = M_sr; // начальное условие для вычислений уравнения Кеплера
      for (;;)   // находим эксцентрическую аномалию через уравнение Кеплера
         {
            E1 = M_sr + (GP[x].E_a * sin(E)); // Kepler_equation(x);
            if(fabs(E1 - E) > 1e-19)
               E = E1;                       // E это по документу  Ek
            else
            break;
         }
        ////////////////////////коррекция синхроимпульсов
       double F,d_tr, d_tsv;
       F = - 4.442807633E-10;                                  // (-2 * sqrt( gp_earth)) / (c * c);
        d_tr = F * GP[x].E_a * GP[x].square_big_axis * sin (E); // E_a - это бывшая e
        d_tsv = GP[x].af0 + (GP[x].af1*(t - GP[x].t_oc)) + (GP[x].af2 * pow (t -GP[x].t_oc,2)) + d_tr - GP[x].t_gd;
                           ;
      // t = t - ( PR_sv[x] / c );//
        t = t - d_tsv; // вычисляем t с учетом поправки
        t_k = t - GP[x].t_oe ;      // вычисляем промежуток времени от опорной эпохи эфемерид
       if ( t_k > 302400)      // если t больше 302400 секунд, вычтите 604800 секунд из T.
        t_k = t_k - 604800;
        else if ( t_k < -302400) // Если T меньше -302400 секунд, добавьте 604800 секунд к T.
        t_k = t_k + 604800;
 
       M_sr = GP[x].Mo + (n * t_k); // определяем среднюю аномалию
 
        E1 = M_sr;     // начальное условие для вычислений уравнения Кеплера
        for (;;)   // находим эксцентрическую аномалию через уравнение Кеплера
         {
           E1 = M_sr + (GP[x].E_a * sin(E)); // Kepler_equation(x);
           if(fabs(E1 - E) > 1e-19)
             E = E1;                           // E_a - это бывшая e
           else
             break;
         }
    v = (atan((sqrt((1 + GP[x].E_a)/(1 - GP[x].E_a))) * tan(E/2)))*2;
    //v = atan((sqrt(1-pow(GP[x].E_a,2))*sin(E)/(1-GP[x].E_a*cos(E)))/((cos(E)-GP[x].E_a)/(1-GP[x].E_a*cos(E))));
             printf("v : %.12f \n",v);
    nash = v + GP[x].omega; // находим невозмущенный аргумент широты
            
    sigma_u = GP[x].Cus * sin(2 * nash) + GP[x].Cuc * cos(2 * nash);
    sigma_r = GP[x].Crs * sin(2 * nash) + GP[x].Crc * cos(2 * nash);
    sigma_i = GP[x].Cis * sin(2 * nash) + GP[x].Cic * cos(2 * nash);
     
    u_eph = nash + sigma_u;
    r_eph = a_bp * (1 - GP[x].E_a * cos(E)) + sigma_r;
    i_eph = GP[x].io + sigma_i + GP[x].IDOT * t_k;
 
    X1 = r_eph * cos(u_eph);
    Y1 = r_eph * sin(u_eph);
           
    OMEGA_r = GP[x].omega_0 + (GP[x].OMEGADOT - OMEGA_zemli)*t_k - OMEGA_zemli*GP[x].t_oe;
          
    X_wgs84 = X1 * cos(OMEGA_r)- Y1 * cos(i_eph) * sin(OMEGA_r);
    Y_wgs84 = X1 * sin(OMEGA_r)+ Y1 * cos(i_eph) * cos(OMEGA_r);
    Z_wgs84 = Y1 * sin(i_eph);
 }

Вот такие результаты:

X23 : 6582196.1679999996 из файла sp3
X : 6582196.4182939604 результат без поправок
X : 6582196.1691133399 результат с поправками

Y23 : 17468760.1230000000
Y : 17468759.9883824740
Y : 17468760.3942986840

Z23 : 19054289.1160000000
Z : 19054289.1614533770
Z : 19054288.8624326470
Мне хотелось бы понять, - правильно ли я применяю поправки, и до какой точности нужно добиваться совпадения результата с файлом SP3/ желательно не отправлять меня в англоязычный сегмент интернета, а просто ответить.......

 

Добрый день. кажется написал расчет координат по методу наименьших квадратов, - для 4 спутников.
получились вот такие результаты:
координаты спутников: //24/02/19 10:30
SV: 1
X : -13511943.4625868720
Y : 22537777.8291581910
Z : -2579689.0497690514
- * ---- end ---- * -

SV: 7
X : -2361138.8937401767
Y : 21700294.4831329210
Z : 15045482.4396653310
- * ---- end ---- * -

SV: 8
X : -10765179.0955109730
Y : 11387996.6558811390
Z : 21462376.4363597040
- * ---- end ---- * -

SV: 30
X : 7438780.3448691871
Y : 15023634.4812066070
Z : 20578692.6607901340
- * ---- end ---- * -

//Time gps: 37800 псевдодальности:
//SV: 1 Q: 7 Psr: 25109839.64387868
//SV: 7 Q: 4 Psr: 21056670.78789259
//SV: 8 Q: 6 Psr: 21555412.29389584
//SV: 30 Q: 4 Psr: 20686354.37529096

результат:
L : 83.229887854074036 B : 54.286477535936783 H : 107614.049749834040000
Координаты: E : 83,13,47 N : 54,17,11

нет ли возможности проверить правильность моих вычислений?

 

Не рановато координаты определять? Поправки уже все введены?

 

Вот я выше немного как раз об этом и спрашивал, - и код приложил. Все ли поправки я учел, если нет, - что упустил?

 

Добрый день. В рекомендованном исходнике программы на Паскале, есть процедура, вычисляющая азимут и угол возвышения. Не мог кто нибудь написать формулу, по которой здесь угол возвышения считается. Формул в интернете много, - но хотелось бы опознать эту. Вопрос возник оттого, что как мне кажется, я неверно перевожу этот, на первый взгляд простенький алгоритм на Си.

Код:

procedure calcAzEl(Xs,              {satellite ECEF XYZ}
                   Xu  :  vec3;     {user ECEF XYZ}
               var Az,              {azimuth [rad]}
                   El  : real;      {elevation [rad]}
               var stat: boolean);  {calculation succeeded: stat = true}
 
var
   R, p, x, y, z, s: real;
   e: array[1..3,1..3] of real;
   i, k: integer;
   d: vec3;
 
begin
 
x := Xu[1];
y := Xu[2];
z := Xu[3];
p := sqrt(x*x + y*y);
if p = 0.0 then stat := false else stat := true;
 
if stat then
   begin
 
   R := sqrt(x*x + y*y + z*z);
   e[1,1] := - y / p;
   e[1,2] := + x / p;
   e[1,3] := 0.0;
   e[2,1] := - x*z / (p*R);
   e[2,2] := - y*z / (p*R);
   e[2,3] := p / R;
   e[3,1] := x / R;
   e[3,2] := y / R;
   e[3,3] := z / R;
 
   for k := 1 to 3 do
      begin
      d[k] := 0.0;
      for i := 1 to 3 do d[k] := d[k] + (Xs - Xu) * e[k,i]
      end;
 
   s := d[3] / sqrt(d[1]*d[1] + d[2]*d[2] + d[3]*d[3]);
   if s = 1.0 then El := 0.5 * pi else El := arctan(s / sqrt(1.0 - s*s));
 
   if (d[2] = 0.0) and (d[1] > 0.0) then Az := 0.5 * pi else
      if (d[2] = 0.0) and (d[1] < 0.0) then Az := 1.5 * pi else
         begin
         Az := arctan(d[1] / d[2]);
         if d[2] < 0.0 then Az := Az + pi else
            if (d[2] > 0.0) and (d[1] < 0.0) then Az := Az + 2.0 * pi
         end;
 
   end;
 
end; {procedure calcAzEl}

 

Elevation and azimuth computation

З.Ы. То, что малограмотные переводчики с англицкого называют углом возвышения в русскоязычной астрономической литературе называется просто высотой. Высота и азимут – это две координаты в горизонтной (горизонтальной) сферической системе координат.

 

Доброе утро. Докопался до коррекции псевдодальностей, - {correct pseudorange}Pcor[prn] := Praw[prn] + dTclck * c - dTiono * c - dRtrop + Cr
Praw[prn] - псевдодальность полученная от модуля.
dTiono - ионосферная поправка, - которую сейчас считаю,
с - скорость света.
вот где мне взять еще эти три поправки. dTclck , dRtrop, Cr В основном документе я их не увидел.

 

Тропосферную модель для расчета поправки мне подсказали, dTclck вроде бы высчитывается при расчетах методом наименьших квадратов, одновременно с координатами, А вот Cr - так и не понятно, что это и для чего.