JP6027824B2

JP6027824B2 - 衛星測位のための軌道パラメータ生成方法、生成プログラム、並びに、衛星測位システム

Info

Publication number: JP6027824B2
Application number: JP2012197473A
Authority: JP
Inventors: 智哉柴田; 裕昭前田
Original assignee: Lighthouse Technology and Consulting Co Ltd
Current assignee: Lighthouse Technology and Consulting Co Ltd
Priority date: 2012-09-07
Filing date: 2012-09-07
Publication date: 2016-11-16
Anticipated expiration: 2032-09-07
Also published as: JP2014052306A

Description

本発明は、全地球測位システム（ＧＰＳ）に代表される衛星測位システムから送信される衛星測位信号生成技術に関し、より具体的には、衛星測位のための軌道パラメータの生成方法、生成プログラム等に関する。

衛星測位システムは、人工衛星と地上設備とからなる。人工衛星からは、測位信号が送信される。

そのような衛星測位システムの例が、全地球測位システム（ＧＰＳ）である。一般に、ＧＰＳは、それぞれ１５７５．４２ＭＨｚ、１２２７．６ＭＨｚ、および１１７６．４５ＭＨｚを中心とするＬ１、Ｌ２、およびＬ５等と称される複数の周波数を使用して動作する。これらの信号のそれぞれが、それぞれの拡散信号によって変調される。また、衛星測位システムの一事例として、日本において開発されている準天頂衛星システム（ＱｕａｓｉＺｅｎｉｔｈＳａｔｅｌｌｉｔｅＳｙｓｔｅｍ：ＱＺＳＳ）が挙げられる（非特許文献１）。ＱＺＳＳもＧＰＳと同様に、それぞれ１５７５．４２ＭＨｚ、１２２７．６ＭＨｚ、および１１７６．４５ＭＨｚを中心とするＬ１、Ｌ２、およびＬ５などの複数の周波数を使用して動作する方針で開発が進められようとしている。

なお、２０１１年９月３０日には、今後の日本の準天頂生成システム（ＱＺＳＳ）についての閣議決定がなされている（非特許文献２）。非特許文献２には、広義の安全保障に資するために、２０１０年代後半には４機の衛星を保有し、さらに将来は７機体制とすることが記載されている。

一方、利用者側の測位は、複数の人工衛星から送信される測位信号を受信して処理することにより行われる。その際、第一に、受信機の時刻タイミングで生成した拡散符号と受信した信号の拡散符号との位相差を計測することで測位衛星と受信機との間の擬似距離が測定される。第二に、測位信号に重畳されている航法メッセージを解読し、そこから航法メッセージの一部である「軌道パラメータ」を抽出し、抽出された軌道パラメータから人工衛星の位置が計算される。上述の第一の測定と、第二の計算から得られた人工衛星の位置とから、受信機の中で位置計算が行われる。

上述の軌道パラメータは周期的に更新される。米国ＧＰＳの場合は、一般に２時間ごとに更新され、我が国のＱＺＳＳの場合は１時間ごとに更新される。この更新周期は、長ければ長いほど、受信機での再取得処理の手間が軽減することができる。一方で、更新周期が短ければ短いほど、受信機は頻繁に航法メッセージの中にある軌道パラメータを抽出しなければならず、演算処理上の負担が増大する。

ここでいう軌道パラメータは、ＧＰＳの場合、非特許文献３のTable 20-IIによって定義されている下表１に示すＩＯＤＥを除く１６個のパラメータで構成される。

また、時刻及び位置の組である（Ｔ，Ｘ，Ｙ，Ｚ）というデータが軌道のカルテジアン軌道要素表現と呼ばれるのに対して、表１に示すデータのような表現方法をケプラリアン軌道要素表現と呼ばれる。

一般に、ケプラリアン軌道要素とは、時刻Ｔに加えて、Ｍ_O、ｅ、Ａ、Ω、ｉ、ωとを合わせて計７個の軌道要素（パラメータ）のことを指す。一方、ＧＰＳ等の衛星測位システムでは、人工衛星の軌道を精密に表現しなくてはならないから、更に９個の軌道要素（パラメータ）が追加され、計１７個となっている。本明細書においては、これを『拡張型軌道パラメータ』と呼ぶことにする。

図８に、人工衛星の楕円軌道モデルにおける軌道長半径ａと離心率ｅとの関係を示す。図８において、楕円の中心Ｏと地球の位置Ｆ₁との距離は、ａｅで表される。

図９に、人工衛星の楕円軌道モデルにおける平均近点離角Ｍの概念を示す。平均近点離角（mean anomaly）は、軌道運動を行う人工衛星（天体）のある時刻における位置を表すパラメータの１つであり、人工衛星（天体）が近地点を通過してからの時間を軌道周期に対する割合で表される。

図１０に、人工衛星の楕円軌道モデルにおける昇交点赤経Ω、軌道傾斜角ｉ、及び近地点引数ωの概念を示す。昇交点赤経（Right Ascension of Ascending Node）Ωは、惑星を周回する人工衛星（天体）において軌道が赤道面を南側から北側に横切る位置(昇交点)の赤経である。軌道傾斜角（Orbital Inclination）ｉは、惑星を周回する人工衛星（天体）において惑星の赤道面と軌道面がなす角度である。また、近地点引数（Argument of Perigee）ωは、天体の運動する方向に沿った昇交点から近点までの角度である。

一般に、カルテジアン軌道パラメータで表現された衛星位置の時刻暦データから、これを近似するケプラリアン軌道パラメータの算出は、最小二乗法若しくは重み付き最小二乗法によって処理される。以下、一例を挙げる。

まず、ある時刻ｔの時の人工衛星の地球固定座標系の位置ｐをｐ（ｔ）と表すものとする。いま、複数のｐ（ｔ）のセットとしてのカルテジアン軌道パラメータが、
であるとし、ケプラリアン軌道パラメータが、
であるとした時に、以下の関係が成立する。
ただし、ｆは、Ｋを変数とするある関数である。

ここで求めたいのは、Ｋである。そこで、以下の演算処理を繰り返すことで求める。この演算処理方法は、ニュートン法と呼ばれる非線形最小二乗法である。

・・・・・・・

また、ＫがＧＰＳやＱＺＳＳで送信している上述の『拡張型軌道パラメータ』である場合も、従来技術によれば、上述した最小二乗法若しくは重み付き最小二乗法によって、ＰからＫが求められていた。
ここで、表１では
ではなく
としているが、これは単位の違いであり、本発明に置いてはどちらでも構わない。

特開平１１−０６４４８１号公報

宇宙航空研究開発機構："準天頂衛星システムユーザインタフェース仕様書（ＩＳ−ＱＺＳＳ）１．１版"、２００９年７月３１日、インターネット＜URL:http://qzss.jaxa.jp/is-qzss/＞平成２３年９月３０日閣議決定「実用準天頂衛星システム事業の推進の基本的な考え方」、インターネット＜URL: http://www.kantei.go.jp/jp/singi/utyuu/pdf/kakugi_jun.pdf＞ IS-GPS-200＜URL：http://www.gps.gov/technical/icwg/IS-GPS-200F.pdf＞

しかしながら、上述したＱＺＳＳは新しく進化しようとしており、使用される人工衛星についても、従来の準天頂衛星に加えて静止衛星が採用されることとなっている。従来、準天頂衛星の軌道は、概ねｉ（軌道傾斜角）＝４５°、ｅ（離心率）＝０．０９といったものであったので、上述した最小二乗法若しくは重み付き最小二乗法によって、ＰからＫを求めたとしても大きな問題は発生しなかった。この場合の精度を表すグラフを図１１に示す。時間の経過に伴う精度誤差は、＋０．１〜＋１．１ｃｍの間に収まっている様子が分かる。

一方、ＱＺＳＳにも採用が予定されている静止衛星の軌道は、

であり、従来の算出方法で軌道パラメータを求めると、図１２に示すように、最大５ｍもの精度誤差が生じてしまうという課題があった。

この課題について考察するに、仮にｉ（軌道傾斜角）＝０°だとすると、軌道が赤道面を交差しないので、Ω（昇交点赤経）が不定となることが考えられる。また、仮にｉ（軌道傾斜角）がゼロに近い値だとすると、Ω（昇交点赤経）については計算上求まるにしても、数値計算が不安定であるし、Ω_dot（昇交点赤経の時間微分）が２時間（軌道パラメータの更新周期）の間で一定とは限らないという問題点が考えられる。

また、仮に、ｅ（離心率）＝０だとすると、軌道は真円となり、近地点や遠地点を定義できず、Ｍ（平均近点離角）が不定となる。また、仮に、ｅ（離心率）がゼロに近い値だとすると、Ｍ（平均近点離角）は計算上求まるにしても、数値計算が不安定であるし、Ωdotが、２時間という軌道パラメータの更新周期の中で一定とは限らないという問題点が考えられる。

そこで、従来の他の方法を採用することもできるが、以下に述べる通り、上記課題を完全に解決できるものとなっていない。

例えば、ＳＢＡＳ（Satellite Based Augmentation System）で採用されているように、ｘ，ｙ，ｚという（人工衛星の）位置を送信するという方法がある。しかしながら、この方法では、受信側のアルゴリズムの大幅な改変が必要となることに加え、頻繁な情報更新を必要とする。例えば、更新のタイミングをＧＰＳやＱＺＳＳと同様に１時間程度とした場合には、図１３に示すように、最大で４ｍ近い誤差を生じる。

また、別の方法として、特許文献１に開示された特異点回避軌道要素を送信するものもある。

すなわち、特許文献１には、測位衛星から送信された軌道情報を用いて当該測位衛星の軌道位置を計算すると共に、当該測位衛星から電波が伝搬される時間に基づいて当該測位衛星から計測すべき地球上の物体までの距離を計算し、前記軌道位置及び前記距離とを用いて前記物体の位置を計測する測位衛星を用いた測位システムにおいて、前記軌道情報を構成する、軌道位置決定のための軌道要素として、少なくとも軌道傾斜角ベクトルのｘ成分、軌道傾斜角ベクトルのｙ成分、衛星軌道の離心率ベクトルのｘ成分、衛星軌道の離心率ベクトルのｙ成分、軌道傾斜角ベクトルのｘ成分変化率、軌道傾斜角ベクトルのｙ成分変化率、離心率ベクトルのｘ成分変化率、離心率ベクトルのｙ成分変化率からなる特異点回避軌道要素を用いることを特徴とする測位衛星を用いた測位システムが開示されている。

しかしながら、特許文献１に開示された技術をＱＺＳＳに採用するには、受信側のソフトウェア改修が大規模なものとなってしまう。

そこで、本発明は、全地球測位システム（ＧＰＳ）に代表される衛星測位システム等に静止衛星が採用されているような場合に、利用者が利用しやすい形式の静止衛星の軌道パラメータを安定的に供給することを目的とする。

本発明は、少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実施される人工衛星の軌道パラメータ生成方法であって、
ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータが与えられたとき、前記パラメータから拡張型特異点回避軌道パラメータＫ´を構成し、前記拡張特異点回避軌道パラメータＫ´のうちのｉ_dotを固定値に拘束し、残りのパラメータを最小二乗法に基づき算出し、その後、拡張型軌道パラメータＫに変換することを特徴とする。

また、本発明は、少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実施される人工衛星の軌道パラメータ生成方法であって、軌道傾斜角ｉが略０である場合（特異点とみなせる場合）に、ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータｐ（ｔ）が与えられたとき、前記パラメータを地球自転軸周りに回転させ、慣性空間座標系へ座標変更してｐ´（ｔ）を算出し、その後、ｐ´（ｔ）を任意の軸で任意の角度で座標回転処理して、ｐ”（ｔ）を算出し、前記ｐ”（ｔ）から拡張型特異点回避型パラメータＫ´を算出し、その後、拡張型軌道パラメータＫに変換することを特徴とする。

本発明に係る衛星測位のための軌道パラメータ生成方法、生成プログラム、並びに、衛星測位システムにより、利用者が利用しやすい形式の静止衛星の軌道パラメータを安定的に供給することができる。

本発明の一実施形態にかかる衛星測位システムの構成を説明する説明図である。本発明の予備的な実施形態にかかる軌道パラメータ生成方法等を使用した場合の位置精度を説明する説明図である。本発明の一実施形態にかかる軌道パラメータ生成方法等を使用した場合の位置精度を説明する説明図である。（４Ａ）は図３に示した誤差の分散を１０日分のデータとして説明する説明図である。（４Ｂ）は４Ａに示した誤差の分散を特定の方向／軸成分に分解したものを説明する説明図である。（４Ｃ）は４Ａに示した誤差の分散を特定の方向／軸成分に分解したものを説明する説明図である。本発明の一実施形態にかかる軌道パラメータ生成方法等における軌道面座標の回転処理手順を説明する処理フローである。本発明の他の実施形態にかかる軌道パラメータ生成方法等における軌道面座標の回転処理手順を説明する処理フローである。本発明の他の実施形態にかかる軌道パラメータ生成方法等を使用した場合の位置精度を説明する説明図である。人工衛星の楕円軌道モデルにおける軌道長半径ａと離心率ｅとの関係を説明する説明図である。人工衛星の楕円軌道モデルにおける平均近点離角Ｍの概念を説明する説明図である。人工衛星の楕円軌道モデルにおける昇交点赤経Ω、軌道傾斜角ｉ、及び近地点引数ωの概念を説明する説明図である。従来の準天頂衛星の軌道を前提とした場合の軌道パラメータを使用して算出した位置精度を説明する説明図である。ＱＺＳＳにも採用が予定されている静止衛星の軌道を想定した場合に従来手法を用いて算出した位置精度を説明する説明図である。従来のＳＢＡＳ（Satellite Based Augmentation System）を採用し更新タイミングをＧＰＳやＱＺＳＳと同様にした場合の位置精度を説明する説明図である。

以下、本発明に係る衛星測位のための軌道パラメータ生成方法、生成プログラム、並びに、衛星測位システムを実施するための形態について、図面を参照しながら説明する。

［実施例１］
図１に、本発明の一実施形態にかかる衛星測位システムの構成を示す。衛星測位システム１００は、大別して、人工衛星１０１ａ〜１０１ｄと、測距受信局１０２と、地上局１０３と、ＧＰＳ受信機１０４とを含む。地上局１０３は、大別して軌道パラメータ処理部１０３１と送信部１０３２とを含むが、本発明の実施形態において説明される処理や演算は、信号送信を除き、地上局１０３内の図示しない計算機によってソフトウェア処理される。図１において、これらの処理系を代表して軌道パラメータ処理部１０３１とする。

地上局１０３内の計算機は、典型的には、キーボード、マウス、タッチパネル等から構成される入力部と、ディスプレイ等の表示部やスピーカ等の音声出力部で構成される出力部と、コンピュータプログラムやデータ等を記憶するためのハードディスク、ＲＡＭ及び／又はＲＯＭ等で構成される記憶部と、コンピュータプログラムによって様々な数値計算や論理演算を行うＣＰＵによって構成される中央処理部４０３と、チップや電気系統等の制御を行うための制御部と、必要に応じて外部ネットワークにアクセスするためのスロットや光通信を行うためのポート、及び通信インタフェースから構成される通信インタフェース部と、計算機内の各モジュールに電源を供給するための電源部とを含み、これらのモジュールは必要に応じて適宜通信バスや給電線によって接続されている。

ＧＰＳ受信機１０４は、地上の位置を計測するための受信機であり、携帯電話等の端末に組み込まれる場合もあれば、車載される場合もある。基本的な動作は、人工衛星１０１ａ〜１０１ｄから送信される航法メッセージ中の軌道情報に基づいて各人工衛星の軌道位置を計算すると共に、各人工衛星から伝播される時刻情報に基づいて各人工衛星までの距離を計算し、これらの軌道位置と距離とを用いて自身の位置を算出する。

本発明の実施例においては、測距受信局１０２で取得したデータを基に、軌道パラメータ処理部１０３１で将来の各時刻における衛星位置ベクトルの時刻暦データを算出し、この衛星位置ベクトルの時刻暦データから、これを近似する軌道パラメータを算出する。算出された軌道パラメータは送信部１０３２及び人工衛星１０１を経由してＧＰＳ受信機１０４に送信される。ＧＰＳ受信機１０４では、受信した軌道パラメータを基に、任意の時刻における衛星位置を算出する。

この一連の処理の流れのうち、衛星位置ベクトルの時刻暦データから、これを近似する軌道パラメータを算出する処理が本発明の対象とする領域である。

次に、本発明に係る衛星測位のための軌道パラメータ生成方法等において使用される軌道パラメータについて、一般的なパラメータの説明を踏まえて詳述する。

まず、一般的な軌道パラメータは、以下の６つのパラメータから構成される。

次に、

の特異点を回避するために一般的に用いられる軌道パラメータは特異点回避型パラメータと呼ばれ、一例として、以下の６つのパラメータからなる（例えば、Oliver Montenbruck他著, "Satellite Orbits: Models, Methods, Applications", Springer, 2000/8出版を参照）。
但し、ｐ及びｑは、軌道傾斜角ベクトルと呼ばれるものである。

人工衛星は不均一な地球重力などの影響を受けて飛行するので、その軌道運動は複雑である。一方、一般的な軌道パラメータを使用する場合であっても、特異点回避型軌道パラメータを使用する場合であっても、地球重力は質点であり、重力は均一である、一般的に２体問題と呼ばれる力学モデルを前提としている。
したがって、一般的な軌道パラメータを使用する場合であっても特異点回避型軌道パラメータを使用する場合であっても、数時間に渡る時間期間の衛星位置をこれらの軌道パラメータ１セットを用いて、衛星測位システムで要求される精度での人工衛星の位置の算出は困難である。

そこで、種々の補正パラメータを含む『拡張型軌道パラメータ』がＧＰＳでは採用されている（非特許文献３，Table20-II）。本発明の一実施形態における補正パラメータは、ＧＰＳで採用されているものと同じ、以下に示す９つのパラメータからなる。

次に、本発明の一実施形態において、『拡張型特異点回避型パラメータ』として、以下のＫ´を採用する。Ｋ´はＧＰＳで採用されている『拡張型軌道パラメータ』のうち、ｅ（離心率）、Ｍ₀（元期の平均近点離角）、Ω₀（元期の昇交点経度）、ｉ₀（元期の軌道傾斜角）、ω（近地点引数）を表３に示したような『特異点回避型パラメータ』に置き換えたものである。Ｋ´をこのように設定することは、本発明の第１の特徴である。

ここで、『拡張型特異点回避型パラメータ』Ｋ´を用いて任意の各時刻ｔにおける衛星の位置ベクトルＰを求める手順を関数ｆとして表せば、関数ｆは以下のように定義される。

ここで、ｆは、＜入力＞、＜出力＞、＜関数＞の項として以下のように整理される。
＜入力＞
・衛星位置の各時刻：ｔ₁，ｔ₂，・・・，ｔ_n
・拡張型特異点回避型パラメータ：

＜出力＞
・各時刻における衛星位置ベクトル：Ｐ

＜関数＞

ここで、ｐｏｓは、時刻ｔと、
を入力として、以下の手順（特異点回避要素の変換及びＩＳ−ＧＰＳユーザアルゴリズム）に基づき計算することができる。

［特異点回避要素の変換］

［非特許文献３ Table 20-IVに示されるＩＳ−ＧＰＳのユーザアルゴリズム］

（※iterationによりＥ_kを解く）

以上より、
が導かれる。

上記のように関数ｆが定義されれば、各時刻における衛星位置ベクトルを基に［背景技術］で示したような、一般的な非線形最小二乗法を用いて『拡張型特異点回避型パラメータ』Ｋ´を求めることができる。

上記の手順で算出された拡張型特異点回避軌道パラメータ
は、次式により、ＧＰＳで採用されている拡張型軌道パラメータ
に変換される。

以上が本発明の実施例１が対象とする図１の軌道パラメータ処理部１０３１で行われる『拡張型軌道パラメータ』Ｋの算出方法である。

本発明により算出された『拡張型軌道パラメータ』ＫはＧＰＳで採用されているものと同一であるため、図１のＧＰＳ受信機１０４では、非特許文献３ Table 20-IVに示される従来と同じ計算方法により、『拡張型軌道パラメータ』Ｋを基に任意の時刻における衛星位置を算出可能である。

以上により求められた『拡張型軌道パラメータ』Ｋを使用した場合の、図１のＧＰＳ受信機１０４で実施される衛星位置算出における、衛星位置精度の様子を図２に示す。図に示されるように、このままでは、４〜５日経過した頃から６ｍの誤差が確認でき、８〜１０日にかけて、８〜１０ｍの誤差が発生する。

そこで、このような誤差の増大を回避するために、Ｋ´においてｉ_dot＝０に設定し、上述と同様の最小二乗法でＫ´を算出する。本発明の一実施形態において、Ｋ´においてｉ_dot＝０に強制設定することが本発明の第２の特徴となっている。ｉ_dot＝０に強制設定されたＫ´は以下のように表される。

ｉ_dot＝０に強制設定されたＫ´を採用した場合の位置精度の様子（３Ｄ誤差）を図３に示す。図に示されるように、誤差は１年を通して２０ｃｍ程度に抑えられており、極めて良好な結果が得られた。

ここで、ｉ_dot＝０に強制設定することの現実的な意義は、ＧＰＳの航法メッセージフォーマットから導かれる。例えば、航法メッセージのビット長が１５ビットであるときには、２の補数表現で、ＬＳＢが２^-44となるので、その最大の範囲は-9.3＊１０^-10〜9.3＊１０^-10となる。そのために、その約１０倍の範囲である-１０^-8〜１０^-8の範囲でｉ_dotを設定することが望ましい。
なお、本実施例では、発明の理解のために、ｉ_dot＝０に強制設定することをターゲットとしたが、当業者から見て等価とみなせる範囲において、ｉ_dot以外の他のパラメータを０以外の他の固定値に強制設定することとしても良いことは言うまでもない。

＜考察＞
図３に示した誤差の分散を短期間分取り出し、特定の方向成分に分解したものを図４（Ａ）、（Ｂ）、及び（Ｃ）に示す。（Ａ）は、図３と同様の３Ｄ誤差の１０日分のデータである。（Ｂ）は、（Ａ）における衛星測位のユーザの測位精度に密接に関わる地球方向（ラジアル方向）であり、（Ｃ）は、（Ａ）における地球固定座標系のＺ軸成分であり静止衛星の場合、クロストラック成分と一致する。図４を参照すると、まず、地球固定座標系のＺ軸成分の誤差が相対的に大きく、図３とほぼ同程度の２０ｃｍ程度の誤差に達するものであり、影響力の大きい誤差要因であることがわかる。

一方、地球方向の誤差は、相対的に極めて小さく１ｃｍ以下であることがわかる。これを衛星測位のユーザ誤差に換算すると、１ｃｍの誤差に加えて、静止衛星を仰角０°で望むような地域においては、ほぼ６分の１程度の誤差となることを意味する。つまり、１ｃｍ＋２０ｃｍ／６＝４ｃｍの誤差である。これは、上述した準天頂衛星の誤差である１ｃｍ程度よりも相対的に大きいものの、ユーザにとっての誤差としては十分許容範囲にあると言える。

＜小括＞
以上説明した本発明の一実施形態によれば、本発明の本質的特徴として次の３つの構成要素を挙げることができる。
（１）新しい『拡張型特異点回避型パラメータ』Ｋ´を導入する。
（２）上記（１）のＫ´において、ｉ_dotを０として設定する。
（３）上記（２）のＫ´を変換処理してＧＰＳで採用されているものと同一の『拡張型軌道パラメータ』Ｋを算出する。

［実施例２］
実施例２では、ｉ＝０に近いという軌道傾斜角の特異点に関する課題を解決するために、事前の座標回転処理を導入する。

まず、地球固定座標系で表現された時刻ｔにおける衛星位置ベクトルｐ（ｔ）を次式により地球自転軸周りに回転させて慣性空間座標系へ座標変更し、それをｐ´（ｔ）とする。その慣性空間座標系とは、ｔ＝ｔ₀の時に自転を止めた慣性空間である。

但し、

：Ｚ軸（地球自転軸）周りの回転行列
である。

その後、慣性空間座標系で表現した時刻tにおける衛星位置ベクトルｐ´(t)を任意の軸で、任意の角度で座標回転させる。ここでは、一実施形態としてとして、Ｘ軸周りにΘだけ回転させるものとする。Θの大きさは、ｉ＝０という特異点を回避するために十分な大きさとすれば良い。ｐ´(t)をＸ軸周りにΘだけ回転させたものをｐ” (t)とすると、
ここで、
：Ｘ軸周りの回転行列
と表すことができる。

次に、このように回転させた時刻ｔにおける衛星位置ベクトルｐ”（ｔ）からなる衛星位置ベクトルの時刻暦データＰ”を使って、Ｐ”が関数ｆ´として表された以下の式を解くことにより（必要に応じて、上述した非線形最小二乗法を用いる）、『拡張型特異点回避型パラメータ』Ｋ´を算出する。実施例１ではｉ_dotを０としてＫ´を算出したが、実施例２においてはｉ_dotを０とはせずに全てのＫ´の要素を算出する。これにより算出された、ｉに相当する変数には、本来の値である
という値ではなく、Θに応じた一定のオフセット角度がついている。Ｐ”及びＫ´は次の通りである。

ここで、『拡張型特異点回避型パラメータ』Ｋ´からＰ”を求める関数ｆ´は実施例１において定義した関数ｆにおいて、地球自転レートをゼロ、すなわち、
としたものである。

上記の手順で算出された拡張型特異点回避軌道パラメータ
は、実施例１と同様にして、ＧＰＳで採用されている拡張型軌道パラメータ
に変換される。

以上の処理が、実施例２において対象とする図１の軌道パラメータ処理部１０３１で行われる『拡張型軌道パラメータ』Ｋの算出方法である。

実施例２により算出された『拡張型軌道パラメータ』Ｋは、ＧＰＳで採用されているものと同一のパラメータ形式であるため、図１の送信部１０３２及び人工衛星１０１では従来と同じデータフォーマットで『拡張型軌道パラメータ』ＫをＧＰＳ受信機１０４に配信可能である。

しかしながら、実施例２における『拡張型軌道パラメータ』Ｋを用いて、非特許文献３ Table 20-IVに示される従来と同じ計算方法では、本来求めるべき衛星位置は計算できない。そこで、実施例２ではＧＰＳ受信機１０４で行うべき新たな計算方法が必要であり、これも本発明の実施例２では含まれる。

ＧＰＳ受信機１０４で行うべき新たな計算方法について説明する。

公開されたユーザアルゴリズム（非特許文献３ Table 20-IV）によれば、軌道面座標での位置を求めた後は、図５に示す処理手順のように座標回転をさせることになっている。そして、通常の受信機には、図５に示すアルゴリズムに実行するソフトウェアが組み込まれており、ＧＰＳ衛星の位置を算出する。すなわち、図５において、軌道面座標での位置情報は、ステップＳ５０２でＲ_x（ｉ_k）という回転処理が施され、ステップＳ５０３でＲ_z（Ω₀＋Ω_dot*tk）という回転処理が施され（この時点で、ＥＣＩ座標系での位置が分かる）、ステップＳ５０４でＲ_z（ＥＲＡ）という回転処理が施されて、ＥＣＥＦ座標系での位置が算出される。

しかしながら、実施例２における処理では、上述の前提により、人工衛星をＸ軸回りにΘだけ回転させた場合の位置が算出される。従って、人工衛星それ自体の位置を知るためには、さらにΘだけ戻す処理を実施する必要がある。なお、この場合の処理に必要な演算要素は定数行列であるので、ユーザアルゴリズムとしては極めて簡便な改修を要するだけである。図６に改修された演算手順を示す。

図６において、軌道面座標での位置情報は、ステップＳ５０２でＲ_x（ｉ_k）という回転処理が施され、ステップＳ５０３でＲ_z（Ω₀＋Ω_dot*tk）という回転処理が施され（この時点で、Ｘ軸にΘｉ回転させたＥＣＩ座標系での位置が分かる）、続いてステップＳ６０４でＲ_x（Θｉ）という回転処理が実施され、ステップＳ６０５でＲ_z（ＥＲＡ）という回転処理が施されて、ＥＣＥＦ座標系での位置が算出される。

上記、『拡張型軌道パラメータ』Ｋを用いて、任意時刻における地球固定座標系での衛星位置ベクトルを算出する手順をより具体的に数式で示す。

＜入力＞
・衛星位置の時刻：ｔ
・拡張型特異点回避型パラメータ：

＜出力＞
・時刻ｔにおける地球固定座標系での衛星位置ベクトル：

＜処理＞

（※iterationによりＥ_kを解く）

より、軌道面座標での衛星位置ベクトル
を次式とする。

軌道面座標での衛星位置ベクトル
を次式により座標変換して慣性空間座標系での衛星位置ベクトル
を求める。

ここで、各軸に沿ってΘ回転させるための定数行列は、以下に示す通りである。

慣性空間座標系での衛星位置ベクトル
を次式により座標変換して最終的に地球固定座標系での衛星位置
を求める。

いま、仮に、Θを３０°とした場合の３Ｄ位置精度を図７に示す。図７に示されるように、誤差は最大でも２ｃｍ程度となり、実施例１で示した手法（最大で２０ｃｍ程度。図３参照）よりも更に良好な結果が得られた。

実施例１では、拡張型特異点回避型パラメータＫ´の1つを固定した為に、３Ｄ誤差で最大２０ｃｍ程度の誤差が残ったが、図１のＧＰＳ受信機１０４で行われる、衛星位置計算方法は従来と同一で構わないことを特徴とする。一方、実施例２では図１のＧＰＳ受信機１０４で行われる、衛星位置計算方法は従来と若干異なるが、３Ｄ誤差で最大２ｃｍ程度と実施例１よりも更に高精度であることを特徴とする。

［公知の技術等］
本発明に関連して、本明細書と同時に出願されたかその前に出願され、公に自由に入手できるすべての論文および文書の内容は、参照によって本明細書の記載内容として組み込まれる。

［組み合わせ］
本明細書（請求項、実施例、要約、及び図面を含む）に記載された構成要件の全て及び／又は開示された全ての方法又は処理の全てのステップについては、これらの特徴が相互に排他的である組合せを除き、任意の組合せで組み合わせることができる。

［特徴の一例］
本明細書（請求項、実施例、要約、及び図面を含む）に記載された特徴の各々は、明示的に否定されない限り、同一の目的、同等の目的、または類似する目的のために働く代替の特徴に置換することができる。したがって、明示的に否定されない限り、開示された特徴の各々は、包括的な一連の同一又は均等となる特徴の一例にすぎない。

本発明は、上述した実施形態のいずれの具体的構成にも制限されるものではない。本発明は、本明細書（請求項、実施例、要約、及び図面を含む）に記載された全ての新規な特徴又はそれらの組合せ、あるいは記載された全ての新規な方法又は処理のステップ、又はそれらの組合せに拡張することができる。

１００衛星測位システム
１０１ａ〜１０１ｄ人工衛星
１０２測距受信局
１０３地上局
１０４ＧＰＳ受信機

Claims

少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実施される人工衛星の軌道パラメータ生成方法であって、
ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータが与えられたとき、前記パラメータから拡張型特異点回避軌道パラメータＫ´を構成し、前記拡張型特異点回避軌道パラメータＫ´のうちの軌道傾斜角の時間微分ｉ_dotを含む１以上のパラメータを固定値に拘束し、残りのパラメータを最小二乗法に基づき算出し、その後、拡張型軌道パラメータＫに変換する
ことを特徴とする方法。
前記固定値は、−１０^-8〜１０^-8の範囲の値に固定されることを特徴とする請求項１に記載の方法。
少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実施される人工衛星の軌道パラメータ生成方法であって、少なくとも軌道傾斜角ｉが略０である場合（特異点とみなせる場合）に、
ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータｐ（ｔ）が与えられたとき、前記パラメータを地球自転軸周りに回転させ、慣性空間座標系へ座標変更してｐ´（ｔ）を算出し、その後、ｐ´（ｔ）を任意の軸で任意の角度で座標回転処理して、ｐ”（ｔ）を算出し、前記ｐ”（ｔ）から最小二乗法に基づき拡張型特異点回避型パラメータＫ´を算出し、その後、拡張型軌道パラメータＫに変換する
ことを特徴とする方法。
前記拡張型軌道パラメータから軸回転座標における位置を求め、定数行列をかけることで前記人工衛星の位置を算出することを特徴とする請求項３に記載の方法。
前記拡張型軌道パラメータは、米国の衛星測位システムであるＧＰＳの仕様に準拠する受信機に実装されている軌道パラメータと同一であることを特徴とする請求項１〜４のいずれか１項に記載の方法。
請求項１〜５のいずれか１項に記載の方法により人工衛星の軌道パラメータを生成し、前記軌道パラメータを測位信号に重畳して前記人工衛星から地上に送信することを特徴とする衛星測位システム。
少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実行される人工衛星の軌道パラメータ生成のためのコンピュータプログラムであって、前記コンピュータプログラムが前記システムにおいて実行されたとき、
ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータが与えられたとき、前記パラメータから拡張型特異点回避軌道パラメータＫ´を構成し、前記拡張型特異点回避軌道パラメータＫ´のうちの軌道傾斜角の時間微分ｉ_dotを含む１以上のパラメータを固定値に拘束し、残りのパラメータを最小二乗法に基づき算出し、その後、拡張型軌道パラメータＫに変換する
ことを特徴とするプログラム。
前記固定値は、−１０^-8〜１０^-8の範囲の値に固定されることを特徴とする請求項７に記載のプログラム。
少なくとも１以上の人工衛星と、地上局と、受信機とを備えたシステムにおいて実行される人工衛星の軌道パラメータ生成のためのコンピュータプログラムであって、少なくとも軌道傾斜角ｉが略０である場合（特異点とみなせる場合）に、前記コンピュータプログラムが前記システムにおいて実行されたとき、
ある時刻Ｔと時刻Ｔにおける人工衛星位置Ｘ，Ｙ，Ｚをと含むパラメータｐ（ｔ）が与えられたとき、前記パラメータを地球自転軸周りに回転させ、慣性空間座標系へ座標変更してｐ´（ｔ）を算出し、その後、ｐ´（ｔ）を任意の軸で任意の角度で座標回転処理して、ｐ”（ｔ）を算出し、前記ｐ”（ｔ）から最小二乗法に基づき拡張型特異点回避型パラメータＫ´を算出し、その後、拡張型軌道パラメータＫに変換する
ことを特徴とするプログラム。