このページは機械翻訳を使用して翻訳されました。最新版の英語を参照するには、ここをクリックします。
惑星暦を使った海洋航行
この例では、惑星の天体暦と地心慣性から地球中心地球固定 (ECI から ECEF) への変換を使用して、船舶の天体航法を実行する方法を示します。
この例では、Mapping Toolbox™ を使用します。aeroDataPackage 関数を使用して、例のデータもダウンロードする必要があります。
この例では、1947 年にコンティキ号が太平洋を横断した遠征隊が辿ったルートを使用します。トール・ヘイエルダールが率いるこの探検隊は、コロンブス以前の時代にポリネシア諸島に南米から来た人々が住んでいたという説を証明することを目的としていた。遠征は101日間かかり、ペルーのカヤオ港からフランス領ポリネシアのラロイア環礁まで航海しました。
注記:この例では、遠征ルートを大まかに再現します。惑星の暦と ECI から ECEF への変換を単純に示すには、多少の自由が必要です。
積載船トラック
astKonTikiData.mat ファイルをロードします。この例の船舶の軌道、速度、および進路が含まれています。このファイルには、さまざまなトラック ポイントの緯度と経度がそれぞれ変数「lat」と「long」に保存されます。変数には、カヤオ港からラロイア環礁までの 1 日あたり 1 つのトラック ポイントに十分なデータが含まれています。このファイルには、各日の船舶速度(ノット/日単位「V」)と針路(度単位「T」)の値も保存されます。
load astKonTikiData
観察構造を作成する
海上縮尺プロセスは、航海士が船舶の緯度と経度を決定するために実行する一連の手順です。これは、アメリカの実用航海士[1]、航海年鑑[2]、天文年鑑補足[3]に記載されている理論に基づいています。このプロセスでは、六分儀、時計、コンパス、航海図から取得した観測データを使用します。観測された各オブジェクトの切片 (p) と真の方位角 (Z) を返します。この例では、観測データを格納するために観測構造配列 obs を使用します。構造体配列のフィールドは次のとおりです。
h:観察者の目の高さ(メートル)。
IC:六分儀のインデックス補正(度)。
P:ローカル周囲圧力(mb)。
T:現地温度(℃)。
年:観測時の現地年。
月:観測時の現地月。
日:観測時の現地日。
時間:観測時の現地時間。
UTC:観測のための協定世界時。年、月、日、時間、分、秒の 6 つの要素のベクトルとして表されます。
Hs:地平線上の天体の六分儀高度(度)。
オブジェクト:測定に使用される天体(木星、海王星、土星など)。
緯度:観測時の船舶の推定緯度(度)。
経度:観測時の船舶の推定経度(度)。
偏角:天体の赤緯(度単位)。
高度:地球の表面から天体までの距離(km)。
GHA:グリニッジ時間角。グリニッジ子午線に対する天体の角度(度単位)。
簡単にするために、すべての測定はボートの同じ場所、同じ六分儀、同じ周囲温度および圧力で行われるものと仮定します。
obs.h = 4; obs.IC = 0; obs.P = 982; obs.T = 15;
遠征隊は1947年4月28日に出発した。その結果、この日付の観測用の構造を初期化します。
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
6 月と 4 月はどちらも 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));
ヘルパー関数 astHelperLongitudeHour を使用して、ローカル時間を UTC に調整します。この機能は、船舶の推定経度に応じて 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');
米国海軍天文台の TAI-UTC (dAT) の値を使用して、ECI から ECEF への変換行列を計算します。
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';
ECEF の位置を使用して、グリニッジ時間角 (GHA) と偏角を計算します。
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 を使用すると、アメリカ大陸とフランス領ポリネシアの視力低下ソリューションを示すより詳細なグラフを取得できます。
if builtin('license','test','MAP_Toolbox') astHelperVisualization(long,lat,longFix,latFix,'Map') end
船がカヤオからラロイアまで航行するにつれて、経度と緯度の相対誤差が蓄積されます。この誤差は六分儀高度の測定における小さな誤差によって発生します。6 月 9 日には、縮小法によって海王星と木星の真の方位角 (Z) が計算されます。海王星と木星の真の方位角は 180 度近く離れています。この違いにより、相対誤差に小さなピークが生じます。ただし、この誤差は依然として削減方法の誤差範囲内にあります。
参考文献
[1] ボウディッチ、N.アメリカの実用航海士。国家地理空間情報局、2012年。
[2] 英国水路局航海年鑑 2012 商用版。パラダイス出版、2011年。
[3] アーバン、ショーンE.とP.ケネス・ザイデルマン。天文年鑑の補足説明。第3版、University Science Books、2013年。
[4] アメリカ海軍天文台
参考
aeroDataPackage
| planetEphemeris
| juliandate
| deltaUT1
| dcmeci2ecef