Как распарсить RINEX OBS File?

2 лабораторная, формулы, которые нам дали в приложении 2(для нахождения координат приемника), приложение 1 для нахождения навигационных координат спутника...

Для нахождения координат приемника я использовал псевдодальномерный метод
Вот мой метод для нахождения координат приемника, расписал все в комментариях:

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

Код:

public class CalcPosition {
    private static final Integer c = 299792458;
    private static final Double eRot = 7.2921151467e-5;
 
    public static Matrix getPosition(List obs, List nav, LocalDateTime dt) throws NullPointerException {
 
        //нахожу по дате проделанные наблюдения
        DataObs rinexObs = RinexObs.searchTimeObs(obs, dt);
        //получаю данные наблюдения спутников
        List observations = rinexObs != null ? rinexObs.getObservations() : null;
        List pr = new ArrayList();
        List xSv = new ArrayList();
        Matrix xMatr = null;
        if (observations != null) {
            for (Observation observation : observations) {
                //нахожу id спутника по наблюдению
                Integer id = RinexObs.getId(observation);
                //ищу в навигационных данных(nav) дату dt и спутник с id
                SatelliteObject navigation = Nav.getNavTimeID(nav, dt, id);
                //в файле присутствует тип спутника R, но я использую только G, если спутник найден
                if (navigation != null && !observation.getSatelliteTypeId().getType().equals("R")) {
                    //получаю PRN, EPOCH,SV CLK
                    SatellitePRN sPRN = navigation.getSatellitePRN();
                    //из наблюдения получая C1
                    BigDecimal c1 = observation.getC1();
                    //получаю время (оно будет равно dt)
                    LocalDateTime epoch = rinexObs.getTime();
                    Long obsTimeMod = (long) (3600 * (24 + epoch.getHour()) + 60 * epoch.getMinute() + epoch.getSecond());
                    BigDecimal tmp = c1.divide(BigDecimal.valueOf(c),11, RoundingMode.HALF_UP);
                    //time = obsTimeMod - C1/299792458; - по формуле
                    Double time = obsTimeMod - tmp.doubleValue();
                    //здесь по времени time и навигационных данных navigation получаю координаты спутник, Вычисление координат на момент излучения
                    ComponentsSatellite cSat = TaskCoordinates.getCoorSat(navigation, null, time);
                    //получаю время из навигационных данных
                    LocalDateTime t = navigation.getSatellitePRN().getTime();
                    Integer toc = 3600 * t.getHour() + 60 * t.getMinute() + t.getSecond();
                    Double tocTime = time - toc;
                    BigDecimal tmp1 = sPRN.getVal("сдвиг часов").add(sPRN.getVal("скорость ухода").multiply(BigDecimal.valueOf(tocTime)));
                    BigDecimal tmp2 = tmp1.add(sPRN.getVal("ускорение ухода").multiply(BigDecimal.valueOf(Math.pow(tocTime, 2))));
                    BigDecimal tmp3 = tmp2.subtract(navigation.getBroadcastOrbitSix().getTgd());
                    BigDecimal dltT = tmp3.add(BigDecimal.valueOf(cSat.getDltTr()));
 
                    //Коррекция псевдодальности, компенсирующая уход часов НКА, т.е приводящая измерения к диной
                    //шкале времени - шкале времени GPS
                    tmp1 = dltT.multiply(BigDecimal.valueOf(c));
                    tmp2 = c1.add(tmp1);
                    //заношу это в лист для сохранения
                    pr.add(tmp2);
                    //Коррекция координат НКА, учитывающая вращение Земли за время от момента излучения, на
                    //который посчитаны координаты НКА, до Момента приём сигнала
                    BigDecimal deltaT = tmp2.divide(BigDecimal.valueOf(c), RoundingMode.HALF_UP);
                    BigDecimal aoR = deltaT.negate().multiply(BigDecimal.valueOf(eRot));
                    //инициалиизруем и обнуляем матрицу mor
                    double[][] mor = getZeroMatr(3, 3);
                    mor[0][0] = Math.cos(aoR.doubleValue());
                    mor[1][1] = mor[0][0];
                    mor[0][1] = Math.sin(aoR.doubleValue());
                    mor[1][0] = -mor[0][1];
                    mor[2][2] = 1.0;
                    Matrix matrix = new Matrix(mor);
                    Matrix transMor = matrix.transpose();
                    double[][] coordSat = getVector(cSat.getSatelliteCoordinates());
                    Matrix coordSatVector = new Matrix(coordSat);
 
                    xSv.add(coordSatVector.times(transMor));
                }
                //Анализ данных на достаточность для обработки
                if (xSv.size() == 4) {
                    break;
                }
            }
            if (xSv.size() < 4)
                return null;
            if (!xSv.isEmpty()) {
                //Первое решение НЗ
                // Формирование начальных значений
                double[][] x = getZeroMatr(4, 1);
                xMatr = new Matrix(x);
                double[][] r = getZeroMatr(4, 1);
                double[][] h = getZeroMatr(4, 4);
                double[][] dPr = getZeroMatr(4, 1);
                double[][] xSt = getZeroMatr(1, 3);
                xSt[0][0] = -2407751;
                xSt[0][1] = -4706536.65;
                xSt[0][2] = 3557571.41;
                for (int i = 0; i < 10; i++) {
                    //Формирование матрицы наблюдений и разницы измеренных и расчётных значений
                    for (int z = 0; z < 4; z++) {
                        r[z][0] = Math.sqrt(Math.pow(getDouble(xSv, z, 0) - xMatr.get(0, 0), 2) + Math.pow(getDouble(xSv, z, 1) - xMatr.get(1, 0), 2) + Math.pow(getDouble(xSv, z, 2) - xMatr.get(2, 0), 2));
                        xSv.get(z).get(0, 0);
                        h[z][0] = (getDouble(xSv, z, 0) - xMatr.get(0, 0)) / r[z][0];
                        h[z][1] = (getDouble(xSv, z, 1) - xMatr.get(1, 0)) / r[z][0];
                        h[z][2] = (getDouble(xSv, z, 2) - xMatr.get(2, 0)) / r[z][0];
                        h[z][3] = 1;
                        dPr[z][0] = pr.get(z).doubleValue() - r[z][0] - xMatr.get(3, 0);
                    }
                    //Вычисление приращения координат
                    //d_X = (H'*H)^-1*H'*d_PR;
                    Matrix matrixH = new Matrix(h);
                    //H'
                    Matrix hTrans = matrixH.transpose();
                    //H'*H
                    Matrix hTransH = hTrans.times(matrixH);
                    //(H'*H)^-1
                    Matrix invHtransH = hTransH.inverse();
                    //(H'*H)^-1*H'
                    Matrix mulInvTransH = invHtransH.times(hTrans);
                    Matrix dPrMatr = new Matrix(dPr);
                    //(H'*H)^-1*H'*d_PR
                    Matrix dX = mulInvTransH.times(dPrMatr);
                    //X = X - dX;
                    xMatr = xMatr.minus(dX);
                }
            }
        } else {
            System.out.println("Не нашло наблюдений спутника");
        }
        return xMatr;
    }
}

 

Какие формулы?
Найдите в сети книгу

Она точно там есть.

Лев Рапопорт до последнего времени сам писал часть firmware для Топконовских приёмников.

Институт, факультет, группа? Кто читает и ведёт практические?

 

Возможно он и хороший преподаватель, но вот за такое "Через центр небесной сферы и Полярную звезду проходит линия называемая осью мира. Точка пересечения оси мира с небесной сферой в направлении Полярной звезды называется Северным полюсом мира, а противоположная – Южным полюсом мира. …Небесный свод вращается вокруг оси мира. При этом Северный полюс мира совмещается с Полярной звездой. "
я бы со спокойной совестью выгнал студента с зачёта. Возможно даже ссаными тряпками гнал. Там столько ляпов по сферической астрономии – закачаяшься. А надо-то было открыть курс общей астрономии Кононовича и Мороза. Во, ещё ляп: "Наклонение орбиты ί – двухгранный угол между экваториальной и орбитальной плоскостями" Во-первых, буква "х" в слове лишняя. Во-вторых, да ты открой школьный курс геометрии, посмотри определение двугранного угла.
"где m = 3,986005×10¹⁴ м³с^-2 - гравитационная постоянная"
Ладно, это всё лирика.
Я вам давал ссылку на GPS Toolkit (www.gpstk.org). В примерах там есть программка PRSolve – читает RINEX (2 и 3 версии), файлы эфемерид (бортовые или точные SP3) и находит навигационное решение. Может это делать по GPS, ГЛОНАСС или по обеим системам. Можно запрещать отдельные спутники, вычислять на разных частотах и т.д.
С вашим RINEX-файлом результат навигационного решения по GPS+ГЛОНАСС на первой частоте (интервал выборки 600 сек, т.е. весовое среднее из 144 значений)
-2407752.1454 -4706538.8351 3557574.2372
из заголовка RINEX
-2407751.0000 -4706536.6500 3557571.4100
на страничке https://earthquake.usgs.gov/monitoring/gps/Southern_California/7odm
-2407751.0427, -4706536.4723, 3557571.6295

 

PRSolve класс очень огромный, сложно там что-то разобрать:(

 

Есть маленько. Хотя, если знаешь немного теорию, то ничего сложного там нет.

 

мне бы только формулы:) Этого было бы достаточно

 

KAI BORRE (раскрыть) KAI BORRE (свернуть)

Gilbert Strang, Kai Borre

Linear algebra, geodesy, and GPS

Ниже работающий код, проверял лет 15 назад.
http://kom.aau.dk/~borre/
http://kom.aau.dk/~borre/life-l99/
В 2012 выходила их книга Algorithms for Global Positioning, но её у меня нет.

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

http://www.insidegnss.com/node/1423
обратите внимание на:
EASY2 computation of a satellite’s position from an ephemeris
EASY3 computation of a receiver’s position from pseudoranges

 

В принципе, практически тоже самое (формулы+код на MathLab в тексте книги) есть в Fundamentals of Global Positioning System Receivers

 

Здравствуйте, а не подскажите алгоритм одномоментного определения координат и отклонения ШВ потребителя по результатам обработки ТНП (текущих навигационных параметров). Где можно это прочитать? гуглил ничего не находит:(

 

Для начала ответьте, почему для получения простейшего решения вам надо четыре спутника, а не три? Почему нельзя получить решение из комбинации 3+1, где 3 - число спутников одной системы (например, GPS), а 1 – число спутников другой (например, ГЛОНАСС)?

 

координаты по идее точнее будут

 

Ещё раз внимательно посмотрите на те уравнения, которые решаете.

 

уравнения решаются методом наименьших квадратов, Находится R и матрица H (4x4)

 

Иногда я бываю назойливее мухи.
Я не спрашиваю как вы решаете уравнения.
Я спрашиваю что за уравнения вы решаете.
Что в этих уравнениях известно и что надо найти.

 

Нам известны координаты навигационных спутников. Нужно найти растояние от приемника до этого самого спутника (R), и координаты этого приемника

 

И…

 

так же известно время

 

Время чего и по каким часам?
Внимательно прочитайте раздел "Общий принцип работы систем спутниковой навигации" из вашей методички и да наступит просветление.

 

аа,

вот в чем проблема, я брал одинаковое время у НС и приемника, а нужно брать разное

 

Оффтоп

Семён Семёнович…