Путь грузового судна

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

load astKonTikiData

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

Процесс морского уменьшения представляет собой серию шагов, которые следует штурман для определения широты и долготы своего судна. Он основан на теории, описанной в «Американском практическом навигаторе» [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;

Инициализация процесса неуплаты для навигации

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

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 в зависимости от предполагаемой долготы. Для этого примера предположим, что все наблюдения берутся одновременно каждый день в 8 вечера по местному времени.

    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 с помощью вспомогательной функции astHelperLegenteHour. Эта функция регулирует время наблюдения UTC в зависимости от расчетной долготы судна.

    obs.UTC = astHelperLongitudeHour(obs);

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

Для каждой из планет astHelperNauticalCalculation функция помощника вычисляет измерение секстанта, которое измерила бы команда в Коне-Тики. Эта функция моделирует фактическое поведение планет, компенсируя при этом локальные условия. Эта функция использует планетарную эфемериду и матрицу преобразования 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')

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

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

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

Ссылки

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

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

[3] Урбан, Шон Э. и П. Кеннет Зайдельманн. Пояснительное дополнение к Астрономическому альманаху. 3-е изд., Университетские научные книги, 2013.

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