天測航法に挑む（前編）

hideです。
　松山空港発・西回り世界一周の途上、リオデジャネイロに寄港して１カ月あまり。日本語サイトも移転を終えて、すっかりきれいになりましたね。tetsuさん、いつもながら本当にありがとうございます。皆様、ご無沙汰を致しました(^^;)。

　実はここしばらく、FlightGearの天測航法（天文航法）を研究していました。toshiさんが先日公開された、天体高度の測定法を基にして、緯度経度を計算する方法を考え、Excelで表を作ってテストを行いました。残念ながら、大小さまざまな誤差が見られ、まだ実用レベルには届きませんが、中間報告させて頂きます。

●天測航法への挑戦：
　六分儀で天体の高度角を測り、緯度経度を算出する天測航法は、航空の世界では1950年代ごろまで、主に渡洋飛行で使われました。現在はINS（慣性航法）やGPSに代替され、船舶でも六分儀の搭載義務がなくなり、天測は過去の技術になりつつあります。しかし航海士の養成課程では実習しますし、海自の哨戒機も、GPS停波などに備えて天測訓練を続けているそうです。
　「太陽や星々を使い、現在地を出す」というのは…なかなかロマンチックな技術です。各種の航海記などを読む際も、天測の知識があると面白くなります。FlightGearは一応、天体の位置を正確に表示していることになっていますので、いつかぜひ、天測航法を研究してみたいと思っていました。

　これまでは画面上で、太陽や星の高度角を測る方法が見つかりませんでしたが、「機体開発」フォーラムの「Nasalを知ろう、使おう!」に先日、「六分儀があると面白いです」と書かせて頂いたところ、toshiさんが月の高度角を実測し、スナップショットを公開（４月25日付）されました。お陰様で、当分は無理かと思われた天測航法は、にわかに現実味を帯びてきました。

　　　　　○

　そこで私としましては、ともかく天体の高度角、方位角、現在時刻（GMT）から、緯度経度を出すExcelシートを作ろうと思いました。
　最初は実際の船舶の、天測航法に関する資料を漁りましたが、後で述べるように、これはFlightGearには、そのままでは転用できないことが分かりました。となると、原理のヒントだけ頂いて、自分で計算法を考えるしかありません。
　しかし私は、数学が非常に（偉大な学問だと思ってますが）苦手です。原理的には、サイン・コサイン、アークタンジェントを使うだけで行ける…と踏んだものの、どこから手を付けたらいいか分からず、自宅の地球儀を見ながら「地軸の傾斜なんて、無ければ簡単なのになぁ」と、ため息をつきました。

　…そうだ。地軸の傾斜角を、無視すればいいのです。
私は、太陽と地球の位置関係を、徹底的に簡素化して、次のような日照モデルを考えました。

・地球は楕円回転体ではなく、真球とする。
・地軸の傾斜角（約23.4度）は、無視する。
・太陽は、赤道の経度ゼロ点上空に、静止している。
　つまり地球の自転も無視する。

…そんなのムチャクチャですって？　いえいえ。グリニッジ時刻の春分の正午は、大体この状態ですから、大丈夫です。まずは単純な特殊例を考え、必要な計算の構造が分かったら、複雑な一般解へ拡張していこう、という要素還元主義的な作戦ですね。

●天体を使って、位置を出す原理：
　ここで「天測航法とは何か」を、おさらいします。理解への突破口は、「位置の線」という概念です。

　いま、あなたが六分儀（または、鉛直方向を示す振り子と、照準器を付けた分度器）で、北極星の高度角を測ったとします。仮に30度だとしますと、あなたは北半球を一周する、北緯30度線上の「どこか」に立っています。こうした、観測者の位置を含む円や線分を、航法では「位置の線」と呼びます。さらに二つの星の高度を測りますと、計３本の位置の線が得られ、その交点（または誤差を見込んで、３本の線に囲まれた三角形の中央）が現在地です。

　昼間は太陽だけを使いますので、一度に１本の位置の線しか得られません。そこで数時間おきに観測をくり返し、その間の推定航行距離だけ、前回得た位置の線を、並行移動して交点を出します。或いは、太陽が南中する瞬間の高度を計って、赤緯（太陽の天球上の位置）などから緯度を算出します。またGMT時間と、南中時の（天文学的な）ローカル時間のずれから、経度を算出します。

　現実の天測では、天体の高度角だけを計ります。先ほどのお話に出て来た、振り子と照準器付きの分度器（大航海時代に使われた、アストロラーベという機械）はもちろん、現代の六分儀にも、方位角を測る装置は付いていません。（コンパスが絡むので、高精度を出すのは無理なのでしょう）
　しかしFlightGearでは、Browse Internal Propertiesの表示を使って、太陽の高度角と方位角、それにGMTの現在時刻が、正確に得られます。これらを使えば１回の観測で、緯度経度の両方を一意に決定できるはずです。

●誤差についての考え方：
　実際の天測航法では、計算式に多くの補正を加えます。例えば天文学的な経年・周年誤差は、毎年更新される、天測暦という膨大な数表で補正します。
　また六分儀は、水平線を基準として、高度角を測る仕組みです。水平線の高さは、観測者の眼高によって変化するため、これを補正します。また気圧と気温、海水温度も計算に加えます（これらは、大気の屈折率に関係があるのでしょうか）。航空機で天測する場合は、風防ガラスの屈折率も問題になるそうです。

　ベテランの航海士が天測をしますと、１nm以内の精度で現在地を測定できるそうです。FlightGearの場合、空港が視界内に入れば十分ですので、10nm前後の誤差まで許容できると思います。
　私は今回、本職の航海士さんが作った、天測計算を自動化するワークシートを入手しました。しかし残念ながら、先に挙げた多くの補正要素が入っているうえ、速度の遅い船舶を前提にした観測・計算法に沿って操作する仕組みで、FlightGearでは当面、あまり使い道がないことが分かりました。

●FlightGear上の観測方法：
　実際にデータを読む方法を、詳しくご紹介します。
　太陽の見える状態でポーズを掛け、視界をHelicopter Viewに切り替えて、HUDを起動します。中央の◇マークを照準器として使い、太陽に合わせます。さらにXキーで太陽をズームアップし、◇マークを太陽の中心に合わせます。
　マウスやタッチパッドの分解能の関係で、◇マークはなかなか、ど真ん中に合ってくれませんが、太陽の視直径は0.5度（角度の30分）ですから、仮に◇が太陽直径の三分の一だけずれると、約10nmの誤差が生じます。私は半透明のUFOでテストしたため、比較的楽に照準できました。慣れたらポーズを掛けず、オートパイロットで飛行しながら測定すると、リアルでしょうね。

　メニューバーの「File」から、Browse Internal Propertiesを開きます。この中にあるsim/current-view/の階層で、以下の数値を読み取ってメモします。
goal-heading-offset-deg：機首を基準とした太陽方位。
　　　　　　　（機首をゼロ度とした、左回り表示です）
goal-pitch-offset-deg：水平線から見た太陽高度。
　　　　　（数値がマイナスですが、符号は無視します）

　またsim/time/の階層から、現在時刻を取得します。
gmt-string：世界時の24時間表示。
　最後にルートの直下にあるorientation/の階層で、機首方位を取得します。
heading-deg：機首方位（真方位）。

（注：Helicopter Viewで観測する場合、天体の高度角は、機体のピッチやロールの影響を受けません。toshiさんの観測で、高度角の測定結果が約6.2度ずれているのは、本来は必要のない/orientation/pitch-deg = 6.25で補正されたからのようですね。このあたりは、確かに分かりにくくて、頭をひねりました）

　…以上で観測はおしまいです。

●緯度経度計算の、基本的な考え方：
　いよいよ緯度経度の計算です。先ほどの「思考モデル」を紙の上に走り書きして、この上に伸縮自在な円形の、位置の線を想像してみました。（マイアルバムに、説明図をアップしました）

　「位置の線」の円の中心点は、赤道上・経度ゼロのアフリカ沖にあります。この位置の線の半径は、太陽の高度角で決まります。別の言い方をしますと…ある地点で太陽の高度角を測ったら、30度だったとします。すると観測者は、北は北緯60度線、南は南緯60度線にそれぞれ接する、円の上のどこかにいるはずです。
（説明図では、青い円に当たります）
　この「どこか」から、太陽が見える方位角を確定すると、現在地が円の上の、ただ１点に決まるわけです。では、ここまでの計算式を作りましょう。

（１）太陽の高度角＝αとする。

（２）太陽の方位角＝βとする。

（３）位置の線上の、任意の観測点をｐとする。
　この、ｐの位置を数式化する方法を考えます。
まず「輪切り」円の中心をｓとし、ｐまでの半径ｒを求めます。ベクトルｓｐの回転角を、ｙ軸を基準にθ度とします。ｐが「輪切り円」（位置の線）上のどこにあろうとも、太陽は常にｓの方向にありますから、太陽の方位角βは常に（θ−180）度で表せます。
　ここでｐからｙ軸までの距離＝ｘと、ｘ軸までの距離＝ｙを考えますと、
　　　ｒ＝sin（90-α）
　　　ｘ＝ｒsinθ
　　　ｙ＝ｒcosθ
という関係になっていることが分かります。

　またｓから、地球の中心ｇまでの長さをＡとしますと、
　　　Ａ＝cos（90-α）
となります。ここでｘを、地球の中心へ伸びる三角形の一辺と考えると、ｐの経度というのは、要するにこの三角形の、地球の中心にある頂角のことだと分かります。同様にｙを一辺とする三角形の頂角が、すなわちｐにおける緯度だと分かります。
　従ってＡを使うと、ｐの緯度経度を、
　　　緯度＝arctan（ｙ/Ａ）
　　　経度＝arctan（ｘ/Ａ）
として算出することができます。

　この計算方法をもう少し発展させ、試行錯誤を重ねた上、ようやくExcelのワークシートを作りました。長文になりましたので、ここまでを前編とします。
（後編は、このあとすぐ!!）