JP6640615B2 - Orbit calculation device and orbit calculation program - Google Patents
Orbit calculation device and orbit calculation program Download PDFInfo
- Publication number
- JP6640615B2 JP6640615B2 JP2016047624A JP2016047624A JP6640615B2 JP 6640615 B2 JP6640615 B2 JP 6640615B2 JP 2016047624 A JP2016047624 A JP 2016047624A JP 2016047624 A JP2016047624 A JP 2016047624A JP 6640615 B2 JP6640615 B2 JP 6640615B2
- Authority
- JP
- Japan
- Prior art keywords
- time
- state quantity
- orbit
- state
- trajectory
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000004364 calculation method Methods 0.000 title claims description 71
- 238000000034 method Methods 0.000 claims description 33
- 239000011159 matrix material Substances 0.000 claims description 29
- 230000007704 transition Effects 0.000 claims description 23
- 230000002123 temporal effect Effects 0.000 claims description 6
- 238000010586 diagram Methods 0.000 description 11
- 230000010354 integration Effects 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004891 communication Methods 0.000 description 3
- 238000013213 extrapolation Methods 0.000 description 3
- 230000001133 acceleration Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 1
- 230000006870 function Effects 0.000 description 1
Images
Landscapes
- Navigation (AREA)
- Complex Calculations (AREA)
Description
本発明は、衛星が定められた周期で同一の位置を通過する基準軌道の近傍を飛行する場合において、観測により取得されたデータに基づいて衛星の軌道を推定・予測・制御する演算処理を行う装置に関する。 The present invention performs an arithmetic process for estimating, predicting, and controlling the orbit of a satellite based on data obtained by observation when the satellite flies in the vicinity of a reference orbit passing the same position at a predetermined cycle. Related to the device.
衛星の利用分野の拡大に伴い、衛星の軌道を推定・予測・制御する性能に対して要求される精度が高くなっている。宇宙から地球の状態をモニタする地球観測衛星においては、観測した地球上の対象位置を特定したり、衛星に作用する加速度を測定したりするために、衛星が飛行した位置を高い精度で推定する必要がある。複数の衛星が定められた間隔で飛行し観測したデータを合成することにより地球上の3次元的な情報を得ることができるため、衛星が定められた位置を飛行するように軌道を予測・制御する技術が求められている。既に市民生活へ浸透している技術であるGPSをはじめとする全地球航法衛星システムにおいては、衛星の電波を受信するユーザが自らの位置を特定するために衛星の軌道情報を利用する必要があり、衛星から配信される軌道情報の精度が測位精度に直接影響する。 With the expansion of satellite application fields, the accuracy required for the performance of estimating, predicting, and controlling the orbit of a satellite is increasing. An earth observation satellite that monitors the state of the earth from space estimates the position where the satellite flew with high accuracy in order to identify the target position on the earth and measure the acceleration acting on the satellite. There is a need. Predict and control the trajectory of a satellite so that it flies at a specified position because it is possible to obtain three-dimensional information on the earth by combining data obtained by multiple satellites flying at specified intervals and observing. There is a need for technology to do this. In global navigation satellite systems such as GPS, which is a technology that has already penetrated civil life, it is necessary for users who receive satellite radio waves to use satellite orbit information to identify their positions. However, the accuracy of the orbit information delivered from the satellite directly affects the positioning accuracy.
衛星の軌道を高い精度で推定する方法として、特許文献1および特許文献2で利用されているバッチ最小二乗推定法が用いられることが多い。この方法を用いる場合には、衛星の軌道力学モデルに基づく数値積分等の計算を繰り返し実施する必要があるため計算機の処理負荷が高く、推定に要する処理時間が長くなるという課題があった。 As a method of estimating the orbit of a satellite with high accuracy, the batch least square estimation method used in Patent Documents 1 and 2 is often used. When this method is used, it is necessary to repeatedly perform calculations such as numerical integration based on the orbital dynamics model of the satellite, so that the processing load on the computer is high, and the processing time required for estimation is long.
この発明の軌道計算装置は、人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する格納部と、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いて、前記人工衛星の状態量の推定値を計算する計算部と
を備えることを特徴とする。
The orbit calculation device according to the present invention provides a reference orbit state which is a state quantity representing a position, a speed, and a dynamic parameter of the artificial satellite in a reference orbit in which the artificial satellite passes through the same position at a predetermined cycle. A storage unit for storing reference trajectory information including a quantity and a reference trajectory state transition matrix that is a state transition matrix representing a temporal transition of a minute change in the reference trajectory state quantity;
A calculating unit that calculates an estimated value of a state quantity of the artificial satellite using the reference orbit information and an observed quantity that is a physical quantity observed for the artificial satellite.
本発明の軌道計算装置では、基準軌道情報をあらかじめ計算し、格納部にデータとして記録し、記録されたデータを利用して、観測量から状態量の推定値を計算する。よって、処理が簡略化され、計算機である軌道計算装置の処理負荷を低減できる。 In the trajectory calculation device of the present invention, the reference trajectory information is calculated in advance, recorded as data in the storage unit, and an estimated value of the state quantity is calculated from the observed data using the recorded data. Therefore, the processing is simplified, and the processing load on the trajectory calculation device, which is a computer, can be reduced.
実施の形態1.
はじめに、一般的なバッチ最小二乗推定法について説明する。
Embodiment 1 FIG.
First, a general batch least squares estimation method will be described.
時刻tにおける衛星の状態量をx(t)と表し、次の式1のように定義する。
状態量x(t)は時刻t0における状態量x0=x(t0)から軌道力学モデルを用いて計算することができる。
状態量x(t)に対する観測モデルをy(t)とする。 The state quantity of the satellite at time t is represented by x (t), and is defined as in the following Expression 1.
The state quantity x (t) can be calculated from the state quantity x 0 = x (t 0 ) at time t 0 using an orbital dynamic model.
Let y (t) be the observation model for the state quantity x (t).
いま、時刻tk(k=1,2,...,n)において観測量z(t)が得られたとき、時刻t0における状態量の推定値はx0に次の式2の変化量Δx0を加算することで得られる。
式3のJk(k=1,2,...,n)はヤコビアンであり、次の式4により計算される。
式4のTijは状態遷移行列と呼ばれ、状態量の時刻tiにおける微小変化δxiに対して時刻tjに生じる微小変化δxjの関係δxj=Tijδxiを表す行列である。
ここで、T0kとT0(k+1)には次の式5の関係がある。
Here, T 0k and T 0 (k + 1) have the relationship of the following Expression 5.
以上のようにして、時刻t0における状態量の推定値が(x0+Δx0)として得られる。その推定値を用いて軌道力学モデルに基づく数値積分等により時刻tにおける状態量x(t)を推定することができる。 As described above, the estimated value of the state quantity at time t 0 is obtained as (x 0 + Δx 0). The state quantity x (t) at time t can be estimated by numerical integration or the like based on the orbital dynamic model using the estimated value.
状態量の推定値を高い精度で得るためには、推定された状態量x(t)を用いて同様の計算を繰り返し、Δx0の各成分の絶対値があらかじめ定めた閾値より小さくなると計算が収束したと判定し、推定処理を終了して状態量の推定値を得る。 In order to obtain the estimated value of the state quantity with high accuracy, the same calculation is repeated using the estimated state quantity x (t), and the calculation is performed when the absolute value of each component of Δx 0 becomes smaller than a predetermined threshold. It is determined that the convergence has occurred, and the estimation process is terminated to obtain an estimated value of the state quantity.
一般的には、時刻t0における状態量x0から時刻tk(k=1,2,...,n)における状態量と状態遷移行列を計算するには軌道力学モデルに基づく数値積分等の方法が用いられ、計算機の処理負荷が高くなることがこの方式の課題となっている。 Generally, in order to calculate a state quantity and a state transition matrix at a time t k (k = 1, 2,..., N) from a state quantity x 0 at a time t 0 , numerical integration based on an orbital dynamic model or the like is performed. The method is used, and the processing load on the computer increases, which is a problem of this method.
図1は、一般的なバッチ最小二乗推定法を示すフローチャートである。次に、図1を用いて一般的なバッチ最小二乗推定法の処理フローについて説明する。 FIG. 1 is a flowchart showing a general batch least squares estimation method. Next, a processing flow of a general batch least squares estimation method will be described with reference to FIG.
ステップS11において、推定処理の開始時に、時刻t0における状態量(初期値)を与える。この時点では状態量の変化量は零であり、ステップS12において、状態量(初期値)を状態量(推定値)とする。 In step S11, providing at the beginning of the estimation process, the state quantity at time t 0 (the initial value). At this time, the change amount of the state quantity is zero, and in step S12, the state quantity (initial value) is set as the state quantity (estimated value).
ステップS12の状態量(推定値)から、軌道力学モデルに基づいた数値積分等を行い、ステップS13で観測モデルyyを計算し、ステップS14でヤコビアンJを計算する。 Numerical integration or the like based on the orbital dynamics model is performed from the state quantities (estimated values) in step S12, the observation model yy is calculated in step S13, and the Jacobian J is calculated in step S14.
ステップS15において、観測により取得された観測量zzからステップS13で計算された観測モデルyyを差し引く。なお観測量zzとは人工衛星に対して観測された物理量である。
ステップS16において、それに対してヤコビアンJと重み行列Wとを用いて、状態量x0の変化量Δx0を計算する。
In step S15, subtracting the observation model y y calculated from the obtained observed quantity z z in step S13 by observation. Note the observed quantity z z is a physical quantity that is observed for satellites.
In step S16, by using the Jacobian J and the weighting matrix W to it, it calculates the amount of change [Delta] x 0 of the state quantity x 0.
ステップS17において、ステップS16で計算された状態量の変化量Δx0の各成分p、v、dの絶対値と閾値とを比較することにより収束判定を行い、閾値より小さい場合には収束と判定し推定処理を終了する。 In step S17, each component p of variation [Delta] x 0 of the calculated state quantity in step S16, v, performs convergence determination by comparing the absolute value and the threshold value of d, the convergence is smaller than the threshold value determination Then, the estimation processing ends.
ステップS17の収束判定において状態量の変化量Δx0の各成分p、v、dの絶対値が閾値より大きい場合には、ステップS12の状態量(推定値)に状態量の変化量Δx0を加算し、それを改めてステップS12の状態量(推定値)として同様の処理を繰り返す。 Each component p of variation [Delta] x 0 of the state quantity in the convergence determination step S17, v, when the absolute value of d is larger than the threshold value, the change amount [Delta] x 0 of the state quantity on the state quantity of the step S12 (estimate) The same process is repeated as the state quantity (estimated value) in step S12.
以下、実施の形態1に係る発明の具体的な実施の形態を説明する。 Hereinafter, specific embodiments of the invention according to Embodiment 1 will be described.
図2は、基準軌道を説明する概略図である。人工衛星が周期Pで同一の位置を通過する軌道は基準軌道と称される。図2は、人工衛星が周期Pで同一の位置を通過する基準軌道40を示している。人工衛星30は、周期Pで位置51を通過し、同様に周期Pで位置52を通過し、同様に周期Pで位置53を通過する。同一の位置とは、位置51、位置52、位置53である。
衛星があらかじめ定められた基準軌道の近傍を飛行する場合には、以下のように推定処理を簡略化することができる。図2では人工衛星30の実際の軌道31を破線で示している。基準軌道40は理想化された軌道である。
FIG. 2 is a schematic diagram illustrating a reference trajectory. The trajectory of the satellite passing through the same position in the period P is called a reference trajectory. FIG. 2 shows a
When the satellite flies in the vicinity of a predetermined reference orbit, the estimation process can be simplified as follows. In FIG. 2, the
人工衛星30が基準軌道40を飛行する場合には、状態量は周期Pで同一の値になる。 このため、
時刻tにおける状態量をx’(t)、
観測モデルをy’(t)と表すと、
x’(t+P)=x’(t)、
y’(t+P)=y’(t)が成り立つ。
なおx’は基準軌道状態量である。
従って、
t(k+m)=tk+Pのとき、
基準軌道40におけるYk、Tk(k+1)をそれぞれY’k、T’k(k+1)と表すと、
Y’(k+m)=Y’k、
T’(k+m)(k+m+1)=T’k(k+1)
が成り立つ。
なおT’は基準軌道状態遷移行列である。
When the
The state quantity at time t is x ′ (t),
When the observation model is represented by y ′ (t),
x ′ (t + P) = x ′ (t),
y ′ (t + P) = y ′ (t) holds.
Note that x 'is a reference trajectory state quantity.
Therefore,
When t (k + m) = t k + P,
When Y k and T k (k + 1) in the
Y ′ (k + m) = Y ′ k ,
T ′ (k + m) (k + m + 1) = T ′ k (k + 1)
Holds.
T 'is a reference orbit state transition matrix.
人工衛星30が基準軌道40の近傍を飛行する場合には、
x(t)≒x’(t)
と近似することができる。
また同様に、
y(tk)≒y’(tk)、
Tk(k+1)≒T’k(k+1)
と近似することができる。
When the
x (t) ≒ x '(t)
Can be approximated.
Similarly,
y (t k) ≒ y ' (t k),
T k (k + 1) ≒ T ′ k (k + 1)
Can be approximated.
時刻t0における状態量の推定値x0は、基準軌道40における状態量x’0に次の式6のように、変化量Δx’0を加算することにより計算する。
つまり、x0=x’0+Δx’0 (式6)
で計算することができる。
式6のΔx’0は、以下の式7、式8、式9、式10で計算する。
That is, x 0 = x ′ 0 + Δx ′ 0 (Equation 6)
Can be calculated by
Δx ′ 0 in Equation 6 is calculated by the following
x’(tk)、y’(tk)、Y’k、T’(k−1)kは、あらかじめ定められた基準軌道40から計算できる。このため、その全て又は一部をあらかじめ計算して計算機の記憶装置にデータとして予め記録しておくことが可能である。また、これらの値は状態量の推定値に依存しないため、繰り返し計算による収束処理を行う必要がない。
x '(t k), y ' (t k), Y 'k, T' (k-1) k can be calculated from the
図3は、上記の式6を模式的に示す図である。上記のように、基準軌道40から計算できるデータを計算機の記憶装置に予め記憶しておき、実施の形態1では、このデータを用いて時刻t0における状態量の推定値x0を計算する。図3において、四角の破線で囲んだデータが、記憶装置に予め記憶されているデータである。これらをまとめて基準軌道情報11Aと呼ぶ。
FIG. 3 is a diagram schematically showing Expression 6 described above. As described above, the data can be calculated from the
観測量を取得する時刻tsにおける
x’(ts)、y’(ts)、Y’s、T’(s−1)sは、
あらかじめ計算された時刻tk(k=1,2,...,n)に対して周期Pの整数(ns)倍を加算しts=tk+nsPが成り立つ時刻tkのデータから得ることができる。
X at time t s to obtain the observables '(t s), y' (t s), Y 's, T' (s-1) s is
Time t k (k = 1,2, ... , n) that has been previously calculated adding the integer (n s) times the period P with respect to t s = t k + n s P holds time t k of the data Can be obtained from
もし、
ts=tk+nsP
が成り立つ時刻tkが存在せず
ts=tk+Δtk+nsP
となる場合には、|Δtk|が最も小さくなる時刻tkを含む1つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。
if,
t s = t k + n s P
T s = t k + Δt k + n s P there is no time tk that is true
And when made is, | Δt k | using data from one or more times, including the smallest becomes the time t k, can be calculated by the method of extrapolation or inner挿等.
以上のようにして、記憶装置に記録したデータである基準軌道情報11Aを用いて推定処理を行うことにより、計算機の処理負荷を低減することができる。
As described above, by performing the estimation process using the
以下に、記憶装置に基準軌道40から得られたx’(tk)、y’(tk)、Y’k、T’(k−1)k等を含む基準軌道情報11Aを予め記憶しておき、基準軌道情報11Aのデータ用いて人工衛星30の軌道の推定値を計算する軌道計算装置10を説明する。
Below, x obtained from the
(***構成の説明***)
図4は、軌道計算装置10及び基準軌道情報生成装置20の機能ブロック図を示す。
図5は、軌道計算装置10及び基準軌道情報生成装置20のハードウェア構成を示す。
軌道計算装置10、基準軌道情報生成装置20はいずれもコンピュータである。
(*** Structure explanation ***)
FIG. 4 is a functional block diagram of the
FIG. 5 shows a hardware configuration of the
The
まず図4の軌道計算装置10を説明する。軌道計算装置10は、格納部11、計算部12、受信部13を備えている。格納部11は、基準軌道情報11Aを格納している。基準軌道情報11Aとは上記のように、基準軌道40における人工衛星30の位置p(t)、速度v(t)、力学パラメータp(t)とを表す状態量である基準軌道状態量x’(t)と、基準軌道状態量x’(t)の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列T’ijとを含む情報である。計算部12は、基準軌道情報11Aと、人工衛星30に対して観測された物理量である観測量zzとを用いて、人工衛星30の状態量の推定値x(t0)を計算する。
受信部13は、後述のように、基準軌道情報生成装置20の生成した基準軌道情報11Aを受信する。格納部11は受信部13が受信した基準軌道情報11Aを格納する。これは格納部11に基準軌道情報11Aがない場合の他、現在の基準軌道情報11Aを新しい基準軌道情報11Aに更新する場合である。
First, the
The receiving
次に図4の基準軌道情報生成装置20を説明する。基準軌道情報生成装置20は、格納部21、生成部22、送信部23を備えている。格納部21は、基準軌道情報11Aを生成するための元データであるオリジナル情報21Aを格納している。生成部22は、オリジナル情報21Aを用いて基準軌道情報11Aを生成する。格納部21が格納している基準軌道情報11Aは生成部22が生成したものである。送信部23は基準軌道情報11Aを軌道計算装置10へ送信する。
Next, the reference trajectory
図5に示すように、軌道計算装置10、基準軌道情報生成装置20はコンピュータである。
As shown in FIG. 5, the
以下では軌道計算装置10を主体に説明するが、この説明は同じコンピュータである
基準軌道情報生成装置20にも当てはまる。軌道計算装置10は、プロセッサ81、メモリ82、通信装置83を備える。プロセッサ81は、プログラムを実行する。メモリ82は格納部11(基準軌道情報生成装置20では格納部21)を実現する。メモリ82には、軌道計算装置10の場合は基準軌道情報11Aが記憶されており、基準軌道情報生成装置20の場合はオリジナル情報21A、基準軌道情報11Aが記憶されている。またメモリ82には、図4に示す「計算部12」(基準軌道情報生成装置20では「生成部22」)の機能を実現するプログラムが記憶されている。そして、プロセッサ81がプログラムを実行して、計算部12(基準軌道情報生成装置20では生成部22)の動作を実行する。通信装置83は受信部13(基準軌道情報生成装置20では送信部23)を実現する。プロセッサ81は、プロセッシングを行うIC(Integrated Circuit)である。プロセッサ81は、CPU(Central Processing Unit)、DSP(Digital Signal Processor)等である。図5に示すメモリ82は、RAM(Random Access Memory)、ROM(Read Only Memory)、HDD(Hard Disk Drive)等である。メモリ82には、計算部12(基準軌道情報生成装置20では生成部22)を実現するプログラムの他にOS(Operating System)も記憶されている。そして、OSの少なくとも一部がプロセッサ81により実行される。図5では1つのプロセッサが図示されているが複数のプロセッサを備えていてもよい。計算部12(基準軌道情報生成装置20では生成部22)の処理の結果を示す情報やデータや信号値や変数値が、メモリ82、又は、プロセッサ81内のレジスタ又はキャッシュメモリに記憶される。
Hereinafter, the
計算部12あるいは生成部22の「部」を、「回路」又は「工程」又は「手順」又は「処理」に読み替えてもよい。また、計算部12あるいは生成部22は、ロジックIC(Integrated Circuit)、GA(Gate Array)、ASIC(Application Specific Integrated Circuit)、FPGA(Field−Programmable Gate Array)といった電子回路により実現されてもよい。なお、プロセッサ及び上記の電子回路を総称してプロセッシングサーキットリーともいう。
The “unit” of the
(***動作の説明***)
図6は、計算部12の処理を示すフローチャートである。図6を用いて、軌道計算装置10の計算部12が格納部11に格納されている基準軌道情報11Aを用いて実行する、「簡略化されたバッチ最小二乗推定法」の処理フローについて説明する。
(*** Explanation of operation ***)
FIG. 6 is a flowchart showing the processing of the
ステップS20において、基準軌道情報生成装置20は、軌道計算装置10による推定処理の開始前に、基準軌道40における人工衛星30の状態量x’、状態遷移行列T’等を計算し、軌道計算装置10に送信する。軌道計算装置10は基準軌道情報11Aを格納部11に記録する。なお、軌道計算装置10が基準軌道情報生成装置20を内蔵する構成でも構わない。つまり、軌道計算装置10が基準軌道情報生成装置20でもある構成でもよい。
In step S20, the reference trajectory
ステップS21において、計算部12は、推定処理の開始時に、格納部11に記録された基準軌道40における状態量x’、状態遷移行列T’等のデータ読み出しを行い、基準軌道40の観測モデルy’yを求め(ステップS22)、また基準軌道40のヤコビアンJ’を求める(ステップS23)。
In step S21, at the start of the estimation process, the
ステップS24において、計算部12は、観測により取得された観測量zzから基準軌道40の観測モデルy’yを差し引く。ステップS25において、計算部12は、その(zz−y’y)に対して、基準軌道40のヤコビアンJ’と重み行列W’を用いて、状態量の変化量Δx’0の計算を行う。
In step S24,
ステップS26において、計算された状態量の変化量Δx’0を、格納部11に記録されている基準軌道40における時刻t0の状態量x’0に加算し、それを状態量の推定値x0(図3)として推定処理を終了する。
In step S26, the calculated change amount Δx ′ 0 of the state quantity is added to the state quantity x ′ 0 at the time t 0 in the
図3、図6に示すように、計算部12は、基準軌道情報11Aを用いて、時刻t0における、状態量の変化量Δx’0を計算し、時刻t0における基準軌道状態量x’0に、計算した変化量Δx’0を加算することにより、時刻t0の推定値x0を計算する。
As shown in FIGS. 3 and 6, the
図7は、軌道計算装置10の基準軌道情報11Aの取得経路を示す図である。図7の(a)は、軌道計算装置10が地上側に存在し、軌道計算装置10が地上側に存在する基準軌道情報生成装置20から基準軌道情報11Aを取得する場合を示している。図7の(b)は、軌道計算装置10が人工衛星30に搭載されており、軌道計算装置10は地上側に存在する基準軌道情報生成装置20から基準軌道情報11Aを取得する場合を示している。なお、軌道計算装置10が基準軌道情報生成装置20を内蔵しても良いことは上記で述べたとおりである。
FIG. 7 is a diagram illustrating an acquisition route of the
(***実施の形態1の効果***)
軌道計算装置10では、基準軌道情報11Aをあらかじめ計算し、メモリにデータとして記録し、記録されたデータを利用したバッチ最小二乗推定法により、観測量から状態量の推定処理を行う。よって、処理が簡略化され計算機の処理負荷を低減できる。
(*** Effect of Embodiment 1 ***)
The
また、軌道計算装置10では、図1のステップS17の判定処理が不要であるので、状態量の推定値を迅速に得ることができる。
In addition, the
実施の形態2.
図8は、実施の形態2を説明するための図である。実施の形態2では、実施の形態1のステップS26で計算された時刻t0の状態量の推定値x0を用いて、計算部12が、時刻t0よりも未来の時刻である時刻tcにおける状態量の予測値xcを計算する。推定値が観測量zzから計算されるのに対して、「予測値」とは、観測量zzを用いることなく計算される状態量の予測結果である。なお実施の形態2は計算部12による処理である。
Embodiment 2 FIG.
FIG. 8 is a diagram for explaining the second embodiment. In the second embodiment, by using the estimated value x 0 of the state of time t 0 which is calculated in step S26 in the first embodiment,
図8に示すように、時刻t0における状態量の推定値x0が得られると、それに基づいて時刻t0よりも未来の時刻tcにおける状態量の予測値xcを、次の式11に従って計算部12が計算する。
時刻tcにおけるx’(tc)、T’(k−1)cは、あらかじめ計算された時刻tk(k=1,2,...,n)に対して周期Pの整数(nc)倍を加算しtc=tk+ncPが成り立つ時刻tkのデータから得ることができる。 X at time t c '(t c), T' (k-1) c is previously calculated time t k (k = 1,2, ... , n) of the period P for integer (n c ) can be obtained from the data at time tk where the times are added and t c = t k + n c P holds.
もし、tc=tk+ncPが成り立つ時刻tkが存在せずtc=tk+Δtk+npPとなる場合には、|Δtk|が最も小さくなる時刻tkを含む1つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。 If the a t c = t k + n c P is absent holds time t k t c = t k + Δt k + n p P is, | one containing the smallest time t k | Δt k Alternatively, it can be calculated by a method such as extrapolation or interpolation using data at a plurality of times.
以上のように計算部12は、実施の形態2によれば、格納部11に記録されたデータと観測量とから推定された衛星の状態量x0を用いて、観測量zzの存在する時刻t0に対して未来の時刻tcにおける人工衛星30の状態量xcを簡便に予測することができる。
以上のように、計算部12は、格納部11に格納された基準軌道情報11Aと、計算された状態量の推定値x0とを用いて、新たな観測量を用いることなく計算される人工衛星30の状態量である予測値xcを計算する。予測値xcは、推定値x0の計算に用いられた観測量zzの存在する時刻よりも未来の時刻tcにおける人工衛星30の状態量である。
As described above, the
実施の形態3.
図9は、実施の形態3を説明するための図である。実施の形態2では予測値xcを計算した。実施の形態3では、計算部12は、時刻tcでの予測値xcを予測値xtにしたい場合の、時刻tf(t0<tf<tc)における状態量の変化量Δxfを計算する場合を説明する。実施の形態3も計算部12による処理である。
Embodiment 3 FIG.
FIG. 9 is a diagram for explaining the third embodiment. In the second embodiment, the predicted value xc was calculated. In the third embodiment,
時刻t0における状態量の推定値x0と時刻tcにおける状態量の予測値xcが得られたとき、時刻tcにおける状態量を目標値xtに変更するためには、時刻tf(t0<tf<tc)において次の式12の変化量Δxfを時刻tfの状態量xfに加える。変化量Δxfは中間変化量である。
時刻tfにおけるT’f(L+1)は、あらかじめ計算された時刻tL(L=1,2,...,n)に対して周期Pの整数(nf)倍を加算しtf=tL+nfPが成り立つ時刻tLのデータから得ることができる。 T ′ f (L + 1) at time t f is obtained by adding an integer (n f ) times the period P to the time t L (L = 1, 2,..., N) calculated in advance, and t f = can be obtained from the data of t L + n f P holds time t L.
もしtf=tL+nfPが成り立つ時刻tLが存在せずtf=tL+ΔtL+nfPとなる場合には、|ΔtL|が最も小さくなる時刻tLを含む1つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。 If the time t L at which t f = t L + n f P does not exist and t f = t L + Δt L + n f P is satisfied, one or more including the time t L at which | Δt L | It can be calculated by a method such as extrapolation or interpolation using data at a plurality of times.
実際に人工衛星30が飛行する軌道を制御する際には、時刻tcにおける状態量の目標値として与えられるのは位置ptあるいは位置ptと速度vtであり、一方、時刻tfにおいて状態量に加えることができるのは速度の変化量Δvfである。
In actual control of the orbit the
実施の形態3によれば、格納部11に記録されたデータと、未来の時刻において予測された衛星の状態量を用いて、その時刻における衛星の状態量を簡便に制御することができる。
According to the third embodiment, the state quantity of the satellite at that time can be easily controlled using the data recorded in the
なお、図9の場合は時刻tfにおける状態量の変化、つまり、状態量(p、v、d)の変化を計算したが、状態量のうちp、vを用いて以下のように計算してもよい。 The change of state quantities at the time t f in the case of FIG. 9, that is, the state quantity (p, v, d) have been calculated changes in, p of the state quantity by using the v calculated as follows You may.
時刻tcにおける状態量の目標値として位置ptが与えられたとき、時刻tfにおいて状態量に加える速度の変化量Δvfは次の式13により計算できる。
時刻tcにおける状態量の目標値として位置ptと速度vtが与えられたとき、2回以上の時刻で状態量に速度の変化量を加える必要がある。時刻tfj(j=1,2,...,N)において状態量に加える速度の変化量Δvfjは次の式15により計算できる。
時刻tfjにおけるT’fjcについても、時刻tfにおけるT’fcと同様の方法で計算することができる。 T ′ fjc at time t fj can be calculated in the same manner as T ′ fc at time t f .
以上のように、計算部12は、基準軌道情報11Aと、状態量の予測値xcとを用いて、中間変化量Δxfを計算する。中間変化量Δxfは、観測量zzの存在する時刻t0よりも未来の時刻tcにおける人工衛星30の予測値xcを異なる予測値xtに変更するための状態量の変化量であり、未来の時刻tcと観測量zzの存在する時刻t0との間の時刻tfにおける状態量の変化量である。
As described above, the
10 軌道計算装置、11 格納部、12 計算部、13 受信部、20 基準軌道情報生成装置、21 格納部、22 生成部、23 送信部、30 人工衛星、40 基準軌道、81 プロセッサ、82 メモリ、83 通信装置。
Claims (7)
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いることにより、時刻t0における状態量の変化量Δx’0を収束判定をすることなく計算し、時刻t0における基準軌道状態量x’0に、計算した前記変化量Δx’0を加算することにより、時刻t0における前記人工衛星の状態量の推定値である推定値x0を計算する計算部と
を備える軌道計算装置。 A reference orbital state quantity, which is a state quantity representing a position, a speed, and a dynamic parameter of the artificial satellite in a reference orbit, which is an orbit through which the artificial satellite passes the same position at a predetermined cycle, and the reference orbital state quantity A storage unit that stores reference trajectory information including a reference trajectory state transition matrix that is a state transition matrix representing a temporal transition of a minute change of
By using the reference orbit information and the observed quantity that is a physical quantity observed for the artificial satellite, the change amount Δx ′ 0 of the state quantity at time t 0 is calculated without convergence determination, and the time t '0, the amount of change Δx calculated' reference trajectory state quantity x at 0 by adding 0, and calculation unit for calculating the estimated value x 0 is an estimate of the state of the satellite at time t 0 Orbit calculation device provided with
前記格納部に格納された前記基準軌道情報と、計算された前記状態量の前記推定値とを用いて、前記推定値の計算に用いられた前記観測量の存在する時刻よりも未来の時刻における前記人工衛星の状態量であり、新たな観測量を用いることなく計算される前記人工衛星の状態量である予測値を計算する請求項1に記載の軌道計算装置。 The calculation unit includes:
Using the reference trajectory information stored in the storage unit and the calculated estimated value of the state quantity, at a time later than the time when the observed quantity used for calculating the estimated value exists. The orbit calculation device according to claim 1, wherein the orbit calculation unit calculates a predicted value which is a state quantity of the artificial satellite and is a state quantity of the artificial satellite calculated without using a new observation quantity.
前記基準軌道情報と、前記状態量の前記予測値とを用いて、前記観測量の存在する前記時刻よりも前記未来の時刻における前記人工衛星の前記予測値を異なる予測値に変更するための状態量の変化量であって、前記未来の時刻と前記観測量の存在する前記時刻との間の時刻における状態量の変化量である中間変化量を計算する請求項2に記載の軌道計算装置。 The calculation unit includes:
A state for changing the predicted value of the artificial satellite at a future time relative to the time at which the observed amount is present to a different predicted value using the reference orbit information and the predicted value of the state quantity; The trajectory calculation device according to claim 2, wherein an intermediate change amount that is a change amount of the amount and that is a change amount of the state amount at a time between the future time and the time at which the observation amount exists is calculated.
前記基準軌道情報を生成する基準軌道情報生成装置から前記基準軌道情報を受信する受信部を備え、
前記格納部は、
前記受信部が受信した前記基準軌道情報を格納する請求項1から請求項3のいずれか一項に記載の軌道計算装置。 The trajectory calculation device,
A receiving unit that receives the reference trajectory information from a reference trajectory information generation device that generates the reference trajectory information,
The storage unit,
The trajectory calculation device according to claim 1, wherein the reference trajectory information received by the reception unit is stored.
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する処理、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いることにより、時刻t0における状態量の変化量Δx’0を収束判定をすることなく計算し、時刻t0における基準軌道状態量x’0に、計算した前記変化量Δx’0を加算することにより、時刻t0における前記人工衛星の状態量の推定値である推定値x0を計算する処理
を実行させるための軌道計算プログラム。 On the computer,
A reference orbital state quantity, which is a state quantity representing a position, a speed, and a dynamic parameter of the artificial satellite in a reference orbit, which is an orbit through which the artificial satellite passes the same position at a predetermined cycle, and the reference orbital state quantity A process of storing reference trajectory information including a reference trajectory state transition matrix which is a state transition matrix representing a temporal transition of a minute change of
By using the reference orbit information and the observed quantity that is a physical quantity observed for the artificial satellite, the change amount Δx ′ 0 of the state quantity at time t 0 is calculated without convergence determination, and the time t '0, the amount of change Δx calculated' reference trajectory state quantity x at 0 by adding 0, executes the process of calculating the estimated value x 0 is an estimate of the state of the satellite at time t 0 Orbit calculation program to make it.
前記基準軌道情報を用いることにより基準軌道のヤコビアンを求め、前記人工衛星に対して観測された物理量である観測量とヤコビアンと重み行列とを用いて時刻t0における状態量の変化量Δx’0を計算し、時刻t0における基準軌道状態量x’0に、計算した前記変化量Δx’0を加算することにより、時刻t0における前記人工衛星の状態量の推定値である推定値x0を計算する計算部と
を備える軌道計算装置。 A reference orbital state quantity, which is a state quantity representing a position, a speed, and a dynamic parameter of the artificial satellite in a reference orbit, which is an orbit through which the artificial satellite passes the same position at a predetermined cycle, and the reference orbital state quantity A storage unit that stores reference trajectory information including a reference trajectory state transition matrix that is a state transition matrix that represents a temporal transition of a minute change of the reference trajectory, and a weight matrix.
The Jacobian of the reference orbit is obtained by using the reference orbit information, and the change amount Δx ′ 0 of the state quantity at time t 0 is obtained using the observed quantity, the Jacobian, and the weight matrix, which are physical quantities observed with respect to the artificial satellite. was calculated, '0, the amount of change Δx calculated' reference trajectory state quantity x at time t 0 by adding the 0, estimated value x 0 is an estimate of the state of the satellite at time t 0 A trajectory calculation device comprising: a calculation unit that calculates
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列と、重み行列とを含む基準軌道情報を格納する処理、
前記基準軌道情報を用いることにより基準軌道のヤコビアンを求め、前記人工衛星に対して観測された物理量である観測量とヤコビアンと重み行列とを用いて時刻t0における状態量の変化量Δx’0を計算し、時刻t0における基準軌道状態量x’0に、計算した前記変化量Δx’0を加算することにより、時刻t0における前記人工衛星の状態量の推定値である推定値x0を計算する処理
を実行させるための軌道計算プログラム。 On the computer,
A reference orbital state quantity, which is a state quantity representing a position, a velocity, and a dynamic parameter of the artificial satellite in a reference orbit, which is an orbit through which the artificial satellite passes the same position at a predetermined period, and the reference orbital state quantity A process of storing reference trajectory information including a reference trajectory state transition matrix, which is a state transition matrix representing a temporal transition of a minute change of, and a weight matrix,
The Jacobian of the reference orbit is obtained by using the reference orbit information, and the change amount Δx ′ 0 of the state quantity at time t 0 is obtained by using the observed quantity, the Jacobian, and the weight matrix, which are physical quantities observed with respect to the artificial satellite. was calculated, '0, the amount of change Δx calculated' reference trajectory state quantity x at time t 0 by adding the 0, estimated value x 0 is an estimate of the state of the satellite at time t 0 A trajectory calculation program for executing the process of calculating the trajectory .
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2016047624A JP6640615B2 (en) | 2016-03-10 | 2016-03-10 | Orbit calculation device and orbit calculation program |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2016047624A JP6640615B2 (en) | 2016-03-10 | 2016-03-10 | Orbit calculation device and orbit calculation program |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2017161420A JP2017161420A (en) | 2017-09-14 |
JP6640615B2 true JP6640615B2 (en) | 2020-02-05 |
Family
ID=59856886
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2016047624A Expired - Fee Related JP6640615B2 (en) | 2016-03-10 | 2016-03-10 | Orbit calculation device and orbit calculation program |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP6640615B2 (en) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR102184662B1 (en) * | 2018-12-17 | 2020-11-30 | 한국항공우주연구원 | A method for predicting satellite events embedded in satellite on-board software |
CN114440886B (en) * | 2021-12-30 | 2023-09-05 | 上海航天控制技术研究所 | High-accuracy track calculation method for large-eccentricity track |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6567712B1 (en) * | 1998-12-02 | 2003-05-20 | Samsung Electronics Co., Ltd. | Method for determining the co-ordinates of a satellite |
JP3550520B2 (en) * | 1999-11-26 | 2004-08-04 | 富士通株式会社 | Trajectory calculation device and trajectory calculation method |
EP1873546A1 (en) * | 2006-06-23 | 2008-01-02 | Nemerix SA | Method and system for ephemeris extension for GNSS applications |
JP2010060489A (en) * | 2008-09-05 | 2010-03-18 | Seiko Epson Corp | Satellite orbital modeling propriety determination method, long term prediction orbital data supply method, and satellite orbital modeling propriety determination device |
EP3124997B1 (en) * | 2014-03-28 | 2020-12-23 | Mitsubishi Electric Corporation | Positioning device |
-
2016
- 2016-03-10 JP JP2016047624A patent/JP6640615B2/en not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
JP2017161420A (en) | 2017-09-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6034615B2 (en) | Vehicle navigation based on satellite positioning data and vehicle sensor data | |
CN107110650B (en) | The method of estimation of controlled navigational state in terms of observability | |
JP5419784B2 (en) | Prediction device, prediction system, computer program, and prediction method | |
Havyarimana et al. | A fusion framework based on sparse Gaussian–Wigner prediction for vehicle localization using GDOP of GPS satellites | |
Ali et al. | Kalman filter tracking | |
US20150153460A1 (en) | Sequential Estimation in a Real-Time Positioning or Navigation System Using Historical States | |
JP2016161570A (en) | Method of obtaining location of device and device | |
JP5692091B2 (en) | Information processing apparatus, information processing method, and computer program | |
JP5151833B2 (en) | Mobile object position estimation system, mobile object position estimation method, and mobile object position estimation program | |
JP2020508527A (en) | Neural episode control | |
JP2020507857A (en) | Agent navigation using visual input | |
KR20210013526A (en) | Apparatus and method for terrain aided navigation using inertial position | |
JP6640615B2 (en) | Orbit calculation device and orbit calculation program | |
WO2016129078A1 (en) | Route selection device and route selection program | |
Hasan et al. | Optimizing of ANFIS for estimating INS error during GPS outages | |
JP5396311B2 (en) | Behavior prediction apparatus, method, and program | |
JP6984739B2 (en) | Search support devices, search support methods, and programs | |
Tharmarasa et al. | Profile-free launch point estimation for ballistic targets using passive sensors | |
Havyarimana et al. | A Hybrid Approach‐Based Sparse Gaussian Kernel Model for Vehicle State Determination during Outage‐Free and Complete‐Outage GPS Periods | |
JP4882544B2 (en) | TRACKING PROCESSING DEVICE, ITS METHOD, AND PROGRAM | |
CN112990595A (en) | Travel time prediction method, travel time prediction device, storage medium and electronic equipment | |
JP6946813B2 (en) | State estimator and program | |
Chen et al. | On SINS/GPS integrated navigation filtering method aided by radial basis function neural network | |
CN113075921A (en) | Local path planning method and device for unmanned equipment | |
JP7491065B2 (en) | State estimation device, state estimation method, and state estimation program |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20180110 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20181024 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20181106 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20181224 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20190611 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20190627 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20191203 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20191226 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 6640615 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
LAPS | Cancellation because of no payment of annual fees |