载重船航迹

加载astKonTikiData。垫文件。在这个例子中，它包含了船的轨迹、速度和航向。该文件分别在变量“lat”和“long”中存储不同轨道点的经纬度。这些变量包含了从卡亚俄港到拉罗亚环礁每天一个航迹点的足够数据。此外，该文件还存储了每天的船舶速度(以节/日“V”为单位)和航向(以度“T”为单位)的值。

负载astKonTikiData

创建一个观察结构

航海还原过程是航海家为确定其船只的经纬度而遵循的一系列步骤。它基于《美国实用导航员》[1]、《航海年鉴》[2]和《天文年鉴》[3]中描述的理论。该过程使用从六分仪、时钟、指南针和航海图获得的观测数据。它返回每个被观测对象的截距(p)和真实方位角(Z)。本例使用观测结构数组obs来包含观测数据。结构数组的字段是:

为简单起见，假设所有测量都是在船上相同的位置，使用相同的六分仪，在相同的环境温度和压力下进行的:

Obs.h = 4;奥林匹克广播服务公司。IC = 0;奥林匹克广播服务公司。P = 982;奥林匹克广播服务公司。T = 15;

探险队于1947年4月28日出发。因此，初始化这个日期的观测结构:

奥林匹克广播服务公司。年份:1947;奥林匹克广播服务公司。月= 4;奥林匹克广播服务公司。天= 28;

初始化导航航迹推算过程

要开始航位推算过程，首先要确定船只位置的初始条件。分别在latFix和longFix变量中存储纬度和经度的固定解决方案的纬度。在本例中，对于第一个固定位置，使用秘鲁卡亚俄的经度和纬度:

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

每日航迹推算

对于本例，假设使用观测数据每天获得一个修正。因此，本例对每个观察使用了“for循环”。变量“m”作为一个计数器，表示自离开港口以来经过的每一天:

为M = 1:size(lat,1)-1

增加日期并对6月和4月的日期进行调整，这两个月都只有30天:

奥林匹克广播服务公司。Day = obs。天+ 1;[obs.month,奥林匹克广播服务公司。day] = astHelperDayCheck(ob .year,ob .month,ob .day);

实际经纬度

从之前装载的航迹点中提取船舶每天的实际位置。下面的例子使用这个值来计算本地时区和所选行星在天空中的位置:

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

地球选择

如果在给定的纬度和经度下，船只可以看到行星，则选择要观测的行星。下面的代码使用预先计算的数据:

如果longActual > -90奥林匹克广播服务公司。对象= {“土星”；“海王星”}；elseiflongActual<=-90 && longActual>-95 obs。对象= {“土星”；“海王星”；“木星”}；elseiflongActual < = -95奥林匹克广播服务公司。对象= {“海王星”；“木星”}；结束

UTC时间计算

根据所假设的经度将本地时间调整为UTC。对于本例，假设所有观测都是在每天同一时间(当地时间晚上8点)进行的。

奥林匹克广播服务公司。小时= 20;

在航迹推算过程中，用当前位置的估计更新观测结构。在这种情况下，位置是使用之前的固定，船只的速度V和航向T来估计的。

奥林匹克广播服务公司。经度= longFix (m) + (1/60) * V (m) *信德(T (m)) / cosd (latFix (m));奥林匹克广播服务公司。纬度= latFix(m)+(1/60)*V(m)*cosd(T(m));

使用helper函数asthelperlongitude dehour将本地时间调整为UTC。该函数根据船只的估计经度调整UTC观测时间。

奥林匹克广播服务公司。UTC = astHelperLongitudeHour(obs);

六分仪高度计算

对于每一个行星，astHelperNauticalCalculation助手函数会计算康提基号上的船员所测量的六分仪测量值。这个函数模拟了行星的实际行为，同时补偿了当地的条件。该函数使用行星星历表和ECI到ECEF转换矩阵。分析不包括行星像差、引力光偏转和光现象的像差。

奥林匹克广播服务公司。Hs = astHelperNauticalCalculation(obs,latActual,longActual);

计算位置

下面的计算取代了航海年历的使用。它们包括行星星历表的使用和ECI到ECEF转换矩阵。

初始化被观测物体的赤纬、格林威治时角(GHA)和高度:

obs.衰落=零(大小(obs.Hs));奥林匹克广播服务公司。GHA = 0 (size(obs.Hs));奥林匹克广播服务公司。height = 0 (size(obs.Hs));

计算测量时间修改后的儒略日:

mjd = mjuliandate(obs.UTC);

计算UT1和UTC之间的差值:

dUT1 = deltaUT1(mjd，“行动”，“没有”）;

使用美国海军天文台的TAI-UTC (dAT)值计算ECI到ECEF的转换矩阵:

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

计算地球时间的儒略日，以近似于地心动力时间:

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

计算每个天体的赤纬、格林尼治时角和高度:

为k = 1:长度(obs.object)

计算每个行星的ECI位置:

posECI = planetEphemeris(jdTT，“地球”, obs.object {k},“405”，“公里”）;

计算每个行星的ECEF位置:

posECEF = TM*posECI';

使用ECEF位置计算格林威治时角(GHA)和赤纬:

obs.GHA(k) = -atan2d(posECEF(2)，posECEF(1));obs. decline (k) = atan2d(posECEF(3)，根号(posECEF(1)^2 +.．.posECEF (2) ^ 2));

使用ECEF to LLA转换函数计算从地球表面到行星中心的距离:

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

结束

行星物体的视力降低

减少观测结构阵列中指定的每个行星的视线:

[p,Z] = elpernauticalreduction (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。经度+ deltaLongFix;latFix(m+1) = obs。纬度+ deltaLatFix;

结束

导航方案可视化

实际履带和视野缩小方案如下图:金宝搏官方网站

astHelperVisualization (lat, longFix latFix,“阴谋”）

您可以使用“绘图工具箱”获取描述美洲大陆和法属波利尼西亚的视力降低解决方案的更详细的图表。金宝搏官方网站

astHelperVisualization (lat, longFix latFix,“地图”）

经纬度的相对误差是随着船从卡亚俄航行到拉罗亚而累积起来的。这个误差是由于六分仪测量高度时的小误差造成的。对于6月9日，还原方法计算出海王星和木星的真实方位角(Z)。海王星和木星的真实方位角接近180度。这种差异导致相对误差的一个小峰值。然而，这个错误仍然在减少方法的错误边界之内。

参考文献

[1]鲍迪奇，N.美国实用航海家。国家地理空间情报局，2012。

[2]英国水道测量局。航海年历2012商业版。天堂出版公司，2011年。

[3]厄本，肖恩·e·p·肯尼斯·塞德尔曼。《天文年历解释补遗》，《大学科学图书》2013年第3版。

[4]美国海军天文台。