Загрузка пути емкости

Загрузите файл astKonTikiData.mat. Он содержит траекторию корабля, скорость и курс для этого примера. Этот файл сохраняет широту и долготу для различных точек дорожки в переменных «lat» и «long» соответственно. Переменные содержат достаточно данных для одной точки дорожки в день от порта Кальяо до атолла Рароя. Кроме того, этот файл также хранит значения скорости каждого дня в узлах в сутки «V» и курсе в град. «T».

load astKonTikiData

Создайте структуру наблюдений

Процесс мореобразования представляет собой ряд шагов, которым следует навигатор для определения широты и долготы своего судна. Он основан на теории, описанной в The American Practical Navigator [1], морском альманахе [2] и пояснительном дополнении к астрономическому альманаху [3]. Процесс использует данные наблюдений, полученные из секстанта, часов, компаса и навигационных графиков. Он возвращает точку пересечения (p) и истинный азимут (Z) для каждого из наблюдаемых объектов. Этот пример использует массив структур наблюдений, obs, чтобы содержать данные наблюдений. Поля для массива структур:

Для простоты предположим, что все измерения проводятся в одном и том же месте в лодке с одним и тем же секстантом при одинаковых температуре и давлении окружающей среды:

obs.h = 4;
obs.IC = 0;
obs.P = 982;
obs.T = 15;

Экспедиция отправилась 28 апреля 1947 года. В результате инициализируйте структуру наблюдения для этой даты:

obs.year = 1947;
obs.month = 4;
obs.day = 28;

Инициализация процесса Dead Reckoning для навигации

Для начала процесса расчетов, задайте начальные условия положения судна. Сохраните широту для фиксированного решения для широты и долготы в переменных latFix и longFix, соответственно. В этом примере для первого местоположения фиксации используйте широту и долготу Кальяо, Перу:

longFix = zeros(size(long));
latFix = zeros(size(lat));
longFix(1) = long(1);
latFix(1) = lat(1);

Ежедневные мертвые расплаты

В данном примере предположим, что исправление получается ежедневно с использованием данных наблюдений. Поэтому этот пример использует цикл for для каждого наблюдения. Переменная «m» действует как счетчик, представляющий каждый истекший день с момента отправления из порта:

for m = 1:size(lat,1)-1

Приращение дня и корректировка дня для месяцев июнь и апрель, оба из которых имеют только 30 дней:

    obs.day = obs.day + 1;
    [obs.month,obs.day] = astHelperDayCheck(obs.year,obs.month,obs.day);

Фактическая широта и долгота

Извлеките фактическое положение судна за каждый день из загруженных ранее точек пути. Пример использует это значение, чтобы вычислить локальный часовой пояс и положение выбранных планет на небе:

    longActual = long(m+1);
    latActual = lat(m+1);

Выбор планеты

Выберите планеты для наблюдения, если они видны судну для заданных широты и долготы. Следующий код использует предварительно вычисленные данные:

    if longActual>-90
        obs.object = {'Saturn';'Neptune'};
    elseif longActual<=-90 && longActual>-95
        obs.object = {'Saturn';'Neptune';'Jupiter'};
    elseif longActual<=-95
        obs.object = {'Neptune';'Jupiter'};
    end

Расчет времени UTC

Настройте местное время по UTC в зависимости от предполагаемой долготы. В данном примере предположим, что все наблюдения берутся одновременно каждый день в 20:00 по местному времени.

    obs.hour = 20;

Для процесса расчета мертвых, обновите структуру наблюдения с оценкой текущего местоположения. В этом случае местоположение оценивается с помощью предыдущего исправления, скорости V судна и курса T.

    obs.longitude = longFix(m)+(1/60)*V(m)*sind(T(m))/cosd(latFix(m));
    obs.latitude = latFix(m)+(1/60)*V(m)*cosd(T(m));

Настройте местное время по UTC с помощью вспомогательной функции astHelperLongitudeHour. Эта функция регулирует время наблюдения UTC в зависимости от расчетной долготы судна.

    obs.UTC = astHelperLongitudeHour(obs);

Расчет высоты секстанта

Для каждой из планет функция помощника astHelperNauticalCalculation вычисляет секстантное измерение, которое экипаж в Kon-Tiki измерил бы. Эта функция моделирует фактическое поведение планет, одновременно компенсируя локальные условия. Эта функция использует планетарную эфемериду и матрицу преобразования ECI в ECEF. Анализ не включает планетарную аберрацию, отклонение гравитационного света и аберрацию световых явлений.

    obs.Hs = astHelperNauticalCalculation(obs,latActual,longActual);

Вычислим положение

Следующие расчеты заменяют использование морского альманаха. Они включают использование планетарных эфемеридов и матрицы преобразования ECI в ECEF.

Инициализируйте склонение, Угол Часа Гринвича (GHA) и высоту для наблюдаемого объекта:

    obs.declination = zeros(size(obs.Hs));
    obs.GHA = zeros(size(obs.Hs));
    obs.altitude = zeros(size(obs.Hs));

Вычислите измененную юлианскую дату для времени измерения:

    mjd = mjuliandate(obs.UTC);

Вычислите различие между UT1 и UTC:

    dUT1 = deltaUT1(mjd,'Action','None');

Вычислите матрицу преобразования ECI в ECEF, используя значения TAI-UTC (dAT) из Военно-морской обсерватории США:

    dAT = 1.4228180;
    TM = dcmeci2ecef('IAU-76/FK5',obs.UTC,dAT,dUT1);

Вычислите юлианскую дату для наземного времени, чтобы аппроксимировать барицентрическое динамическое время:

    jdTT = juliandate(obs.UTC)+(dAT+32.184)/86400;

Вычислите склонение, Угол Часа Гринвича и высоту для каждого из небесных объектов:

    for k =1:length(obs.object)

Вычислите положение ECI для каждой планеты:

        posECI = planetEphemeris(jdTT,'Earth',obs.object{k},'405','km');

Вычислите положение ECEF каждой планеты:

        posECEF = TM*posECI';

Вычислите Угол Часа Гринвича (GHA) и склонение с помощью положения ECEF:

        obs.GHA(k) = -atan2d(posECEF(2),posECEF(1));
        obs.declination(k) = atan2d(posECEF(3),sqrt(posECEF(1)^2 + ...
            posECEF(2)^2));

Вычислите расстояние от поверхности Земли до центра планеты с помощью функции преобразования ECEF в LLA:

        posLLA = ecef2lla(1000*posECEF');
        obs.altitude(k) = posLLA(3);

    end

Сокращение прицела для планетарных объектов

Уменьшите зрение для каждой из планет, заданных в массиве структур наблюдений:

    [p,Z] = astHelperNauticalReduction(obs);

Вычислим шаги к широте и долготе из последнего исправления для текущего исправления используя следующие уравнения. Эти уравнения основаны на морском альманахе.

    Ap = sum(cosd(Z).^2);
    Bp = sum(cosd(Z).*sind(Z));
    Cp = sum(sind(Z).^2);
    Dp = sum(p.*cosd(Z));
    Ep = sum(p.*sind(Z));
    Gp = Ap*Cp-Bp^2;

Вычислим шаги для широты и долготы согласно редукции:

    deltaLongFix = (Ap*Ep-Bp*Dp)/(Gp*cosd(latFix(m)));
    deltaLatFix = (Cp*Dp-Bp*Ep)/Gp;

После вычисления шагов для широты и долготы добавьте их к предполагаемому местоположению, получив исправление на время наблюдения:

    longFix(m+1) = obs.longitude + deltaLongFix;
    latFix(m+1) = obs.latitude + deltaLatFix;

end

Визуализация навигационного решения

Следующий рисунок показывает фактические решения по сокращению дорожки и прицела:

astHelperVisualization(long,lat,longFix,latFix,'Plot')

Можно использовать Mapping Toolbox, чтобы получить более подробный график, изображающий решения по уменьшению прицела с американским континентом и Французской Полинезией.

astHelperVisualization(long,lat,longFix,latFix,'Map')

Относительная погрешность долготы и широты накапливается при плавании судна из Кальяо в Раройю. Эта ошибка связана с небольшими ошибками в измерении высоты секстанта. Для 9 июня метод восстановления вычисляет истинный азимут (Z) для Нептуна и Юпитера. Истинные азимуты Нептуна и Юпитера близки к 180 ° друг от друга. Это различие вызывает небольшой пик в относительную погрешность. Однако эта ошибка все еще находится в пределах контуров ошибки метода сокращения.

Ссылки

[1] Боудич, Н. Американский Практический Навигатор. Национальное управление геопространственной разведки, 2012 год.

[2] Гидрографическое управление Соединенного Королевства. Морской альманах 2012 Commercial Edition. Paradise Publications, Inc. 2011.

[3] Urban, Sean E. and P. Kenneth Seidelmann. Пояснительное дополнение к астрономическому альманаху. 3rd Ed., University Science Books, 2013.

[4] Военно-морская обсерватория США.