Этот пример показывает, как использовать планетарный ephemerides и Землю, В центре Инерционный, чтобы Заземлить Фиксированную Землю В центре (ECI к ECEF) преобразование, чтобы выполнить астронавигацию морского судна.
Этот пример использует Mapping Toolbox™. Необходимо также загрузить данные для примера с помощью aeroDataPackage команды.
Этот пример использует маршрут, сопровождаемый 1 947 экспедициями через Тихий океан Кон-Тики. Экспедиция, во главе с Тором Хейердалем, нацеленным, чтобы доказать теорию, что полинезийские острова были заполнены людьми из Южной Америки в доколумбовы времена. Экспедиция заняла 101 день и приплыла от порта Кальяо, Перу к атоллу Рэроия, Французская Полинезия.
Примечания: Этот пример свободно воссоздает маршрут экспедиции. Требуются некоторые свободы показать планетарный ephemerides и ECI к преобразованию ECEF просто.
Загрузите astKonTikiData.mat файл. Это содержит траекторию поставки, скорость и курс для этого примера. Это хранилища файлов, широта и долгота для различной дорожки указывает в переменных "lat" и "долго", соответственно. Переменные содержат достаточно данных для одной точки дорожки в день от порта Кальяо к атоллу Рэроия. Кроме того, этот файл также хранит значения для скорости судна каждого дня в узлах в день "V" и курс в градусе "T".
load astKonTikiData
Навигационный процесс сокращения является серией шагов, которые навигатор выполняет, чтобы определить широту и долготу его судна. Это основано на теории, описанной в американском Практическом Навигаторе [1], Навигационный Альманах [2], и Объяснительное Дополнение к Астрономическому Альманаху [3]. Процесс использует наблюдательные данные, полученные из секстанта, часов, компаса и навигационных графиков. Это возвращает прерывание (p) и истинный азимут (Z) для каждого из наблюдаемых объектов. Этот пример использует массив структур наблюдения, obs, чтобы содержать наблюдательные данные. Поля для массива структур:
h: Высота уровня глаз наблюдателя, в m.
Ic : индексное исправление Секстанта, в градусе.
P: Локальное окружающее давление, в mb.
T: Локальная температура, в C.
год: Локальный год во время наблюдения.
месяц: Локальный месяц во время наблюдения.
день: Локальный день во время наблюдения.
час: Локальный час во время наблюдения.
UTC: Всемирное координированное время для наблюдения, представленного как шесть векторов элемента с годом, месяцем, днем, часом, минутами и секундами.
Hs: высота Секстанта астрономического объекта выше горизонта, в градусе.
объект: Астрономический объект используется для измерения (т.е. Юпитер, Нептуна, Сатурна, и т.д.).
широта: Предполагаемая широта судна во время наблюдения, в градусе.
долгота: Предполагаемая долгота судна во время наблюдения, в градусе.
наклон: Наклон астрономического объекта, в градусе.
высота: Расстояние от поверхности земли к астрономическому объекту, в км.
GHA: Гринвичский Угол Часа, который является углом в градусах астрономического объекта относительно Гринвичского меридиана.
Для простоты примите, что все измерения проведены в том же местоположении в лодке, с тем же секстантом, при той же температуре окружающей среды и давлении:
obs.h = 4; obs.IC = 0; obs.P = 982; obs.T = 15;
Из экспедиции отбывают 28-го апреля 1947. В результате инициализируйте структуру для наблюдения для этой даты:
obs.year = 1947; obs.month = 4; obs.day = 28;
Чтобы запустить процесс точного расчета, задайте начальные условия для положения судна. Сохраните широту для фиксированного решения для широты и долготы в 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 в зависимости от принятой долготы. В данном примере примите, что все наблюдения берутся одновременно каждый день в 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 функция помощника вычисляет измерение секстанта, которое измерила бы команда в Коне-Тики. Это модели функции фактическое поведение планет, при компенсации локальных условий. Эта функция использует планетарную эфемериду и ECI к матрице преобразования ECEF. Анализ не включает планетарную аберрацию, гравитационное легкое отклонение и аберрацию легких явлений.
obs.Hs = astHelperNauticalCalculation(obs,latActual,longActual);
Следующие вычисления заменяют использование навигационного альманаха. Они включают использование планетарного ephemerides и 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')
Относительная погрешность в долготе и широте накапливается, когда судно приплывает от Кальяо до Raroia. Эта ошибка происходит из-за небольших ошибок в измерении высоты секстанта. На 9-е июня метод сокращения вычисляет истинный азимут (Z) для Нептуна и Юпитер. Истинные азимуты для Нептуна и Юпитер близко к на расстоянии в 180 градусов. Это различие вызывает небольшой пик в относительной погрешности. Эта ошибка, однако, все еще в ошибочных контурах метода сокращения.
[1] Bowditch, N. Американский практический навигатор. Национальная геопространственная спецслужба, 2012.
[2] Гидрографическая служба Соединенного Королевства. Навигационный альманах 2 012 коммерческих выпусков. Paradise Publications, Inc. 2011.
[3] Городской, Шон Э. и П. Кеннет Сейделманн. Объяснительное дополнение к астрономическому альманаху. 3-й Эд., университетские книги по науке, 2013.
[4] Военно-морская обсерватория Соединенных Штатов. https://www.usno.navy.mil.