JP2004150901A - Positioning error simulation system, positioning error evaluation method, program for positioning error simulation, and storage medium for storing the program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a system for performing positioning error simulation with lock delay amount by computing accurate look delay amount in the GPS satellite direction, based on numerical value simulation data expressing almost accurately atmospheric state, and provide a positioning error evaluation method, a program for performing the positioning error simulation, and a storage medium for storing the program. <P>SOLUTION: The positioning error simulation system is equipped with an input means of the numerical value simulation data expressing the atmospheric state, a look delay amount computing means for calculating the look delay amount based on the numerical value simulation data, and a positioning error computing means for obtaining positioning error based on the look delay amount calculated with the computing means. Positioning error is obtained by simulation. <P>COPYRIGHT: (C)2004,JPO

Description

【００１１】
ここで、ＧはＧＰＳ衛星・クエーサーなどの電波源から地上のアンテナまでの直線距離、ｎ（ｓ）は大気屈折によって曲率が生じた伝搬経路Ｌ上の点ｓにおける屈折率である。このとき、Ｓの長さはＬ上の微少部分ｄｓを積分して得られる。右辺第１項が減速の効果であり、実際の電波はＬ上を中性大気の屈折率分だけ遅い伝搬速度で伝搬する。ところが、ＶＬＳＩやＧＰＳの解析の上では同じ経路上を真空中の速度で伝搬したとみなすため、双方の速度差に伝搬時間を乗じた分だけ見かけ上伝搬経路が伸びたことになる。残りの右辺第２項が伝搬経路の曲率の効果となり、当然の事ながらこの効果による遅延は天頂方向で０となる。
【００１２】
ここで、湿潤大気の屈折率の計算は以下のような式で得られる（ ここでは途中の過程を一部省略した） 。
【００１３】
【数２】

【００１４】
ここで、式（２）のＰは全大気圧（ｈＰａ）、Ｐｖ及びＴはそれぞれ水蒸気分圧（ｈＰｗ）、及び絶対温度（Ｋ）である．各項の係数Ｋ_１ 、Ｋ_２ 室内実験より決定される定数であり、研究者によって若干異なる。一般に、式（２）の右辺第１項を“静水圧項”と呼ぶ。この項は、静水圧平衡にある大気成分の密度に比例する寄与であり、これを地上から大気上層まで積分すると静水圧遅延が得られる。式（２）の右辺第２項は“湿潤項”と呼ばれ、これを積分して得られる遅延を湿潤遅延と呼ぶ。
【００１５】
今回の非静力学モデルによるレイトレーシング計算では、まず非静力学モデルの各格子点に与えられた気温、気圧、水蒸気分圧の各値から式（２）を用いて格子点毎の屈折率を計算する。次に、地上のある観測点から任意の仰角と方位で射出した電波が衛星まで到達する経路を、計算した屈折率からＳｎｅｌｌの法則に従って逐次屈折角を計算する。この計算そのものがレイトレーシングである。レイトレーシングで計算された視線方向の遅延量と、従来の大気モデル（これがマッピング関数と呼ばれるもの） とを比較することで、従来のモデルの評価ができる。次に本発明の測位誤差シミュレーションについて説明する。まず従来の大気推定法であるマッピング関数について簡単に説明し、その後でレイトレーシングした結果による測位誤差シミュレーションを説明する。
【００１６】
従来の大気遅延推定、即ち、マッピング関数ＧＰＳやＶＬＳＩの解析では、２つの観測局間での電波の位相差データに含まれる観測局間の相対位置、電離層遅延、大気遅延、観測局ないしは人工衛星内の時計の誤差等の寄与を未知パラメータとして膨大な数の観測方程式を立て、各パラメータを最小二乗法で推定する手法を取る。大気屈折に関していえば、厳密には気温・気圧・水蒸気分圧の時間的・空間的変動が存在するため、視線方向の大気遅延は時間・空間の双方の関数となる。しかし、これらをすべて未知数とするとパラメータが多すぎて観測方程式を解くことが不可能になる。そこで、ある一定の時間内において大気が水平方向に一様な構造を維持すると仮定し、その間は大気遅延が仰角のみに依存するとするモデル化を行えば未知数を大幅に減らすことができる。またモデル化により、デジタル化された観測データを数値計算する上で大気遅延の扱いが容易となる。この思想で導入されたのがマッピング関数と呼ばれるモデルである。マッピング関数を使うと、ある仰角θの全大気遅延ΔＬは下記のように簡単に表すことができる。
【００１７】
【数３】

【００１８】
ここで、ΔＬ^ｚ _ｈ とΔＬ^ｚ _ｗ は、それぞれ天頂方向の静水圧遅延と湿潤遅延である。これらに乗じられているＭ_ｈ （θ）とＭ_ｗ （θ）が、それぞれ静水圧遅延用と湿潤遅延用のマッピング関数である。
【００１９】
後でも触れるが、マッピング関数は、観測方程式を線形化した際に作られる計画行列の中の大気遅延に関する偏微分係数そのものとなる。ただし、ＶＬＳＩでは２地点間でのマッピング関数の差（後述する「一重位相差観測値」に含まれる“大気遅延量差”）、ＧＰＳではさらに異なる衛星について求めた２地点間でのマッピング関数の差どうしで再び差を取った値（後述する「二重位相差観測値」に含まれる“大気遅延量差の差”）となる。
【００２０】
マッピング関数は、基本的には仰角θの正弦の逆数、すなわち（１／ｓｉｎ（θ））となる。ただし、１０度以下の低仰角の観測を行う場合には、マイクロ波が大気中でより長い経路をたどるため、伝搬経路の曲率がきつくなる。その結果、曲率の効果による遅延が指数関数的に増大する。しかも、地球の曲率の影響も無視できなくなる。そこで、単純な「仰角の正弦の逆数」ではなく、下記の式（４）に示すように分母を連分数で拡張したマッピング関数が基本形として考えられている。
【００２１】
【数４】

【００２２】
ここで、各項の係数（ａ、ｂ、ｃ）は定数である。これらの係数は、球対称大気モデルから得られる屈折率プロファイルからレイトレーシングを用いて求めた大気遅延とマッピング関数とを比較し、双方が一致するように最小二乗法を用いて決められる。実際のマッピング関数の開発では、係数を常数としたり、地上気象データの気圧・温度などで変化する変数にしたり、また最後のｂ／（ｓｉｎ θ ＋ｃ）の部分をｔａｎθなど他の関数型に置き換えたりなど、研究者によってより高精度のマッピング関数を求めて様々な改良がなされている。こうして開発されたマッピング関数と、波線追跡法から、それぞれ計算された大気遅延量の差の大小が、マッピング関数の適用範囲、すなわち推定精度を示すことになる。
【００２３】
具体的な手順としては、ＧＰＳやＶＬＳＩ解析で使用される観測方程式に、レイトレーシングした遅延量を入力して観測方程式を解くことにより、大気遅延によって生じる測位誤差シミュレーションが可能になる。
【００２４】
具体的には次のように計算する。まず、ＧＰＳやＶＬＳＩの計算では２地点間の位相差を取ることから、遅延量についても観測点間での差を計算する。観測点Ｒ_１ とＲ_２ に共通の衛星から電波が平行に到達すると仮定すると、観測エポックでｉ番目の２地点間の遅延量差は、
【００２５】
【数５】

【００２６】
ここでΔＬ^Ｓ （Ｒ_１ ）_ｉ とΔＬ^Ｓ （Ｒ_２ ）_ｉ は、それぞれ観測点Ｒ_１ とＲ_２ における任意の方位・仰角での視線遅延量である。添え字のＳとｉはそれぞれ衛星とエポックの番号を示す（ エポックとは観測の区切りである） 。ここで、観測点Ｒ_１ が固定点とすると、遅延量差ΔＬ^Ｓ （Ｒ_２ −Ｒ_１ ）_ｉ は２地点間の遅延量差であり、通常のＧＰＳ解析で言うところの「一重位相差観測値」に含まれる大気遅延量差に相当する。
【００２７】
ここで、多数のエポック毎の遅延量差を観測データとして、この遅延量差に最も適合するベクトルを最小二乗法で計算することができる。すなわち、この関係を線形化すると、
【００２８】
【数６】

【００２９】
である。式（６）で、行列ｅは誤差ベクトル、行列Ｙが遅延量残差から得られる、エポック１番目からｍ番目までの観測量のベクトルとおけて、
【００３０】
【数７】

【００３１】
と表される。また、行列Ａは、各要素の偏微分係数からなる計画行列であり、簡単のため推定するパラメータを観測点位置の誤差ベクトル（後述する見かけの”変位量ベクトル”） と天頂方向の大気遅延量とすると、
【００３２】
【数８】

【００３３】
となる。ここで、
α_ｉ 、β_ｉ およびγ_ｉ （ ｉ＝１、２、・・・・ｍ） などは、遅延量差のｘ、ｙ、ｚの各成分についての方向余弦である。また、Δｍｆ（θ） （Ｒ_２ −Ｒ_１ ）は遅延量についての偏微分係数である（ ここではθを仰角とするマッピング関数の２地点間での差） 。さらに、Ｄが、この方程式であてはめを行った結果得られるベクトルであり、
【００３４】
【数９】

【００３５】
と表現できる。ここで、ｘ、ｙ、ｚはベクトルの各成分であり、ｘ軸の正の方向に北、ｙ軸の正の方向に東、ｚ軸の正の方向に鉛直上向きとする。また、ＺＤは天頂方向の遅延量である。
【００３６】
このベクトルＤは、観測点Ｒ_２ の位置を真の位置からずらすようなベクトルであり、これを見かけの変位量ベクトル、すなわち大気遅延によって生じる測位誤差と見なすことができる。この見かけの変位量ベクトルの計算が測位誤差シミュレーションであり、様々な大気状態の元でこのベクトルの大きさや方向を調べることで測位誤差の振る舞いを調べることが可能になる。
【００３７】
実際に、最小二乗法を用いて遅延量差から見かけの変位量ベクトルＤを求めることを考えると、この場合の観測方程式は、
【００３８】
【数１０】

【００３９】
のように書ける。
【００４０】
ここで得られる行列Ｄは、例えばＶＬＳＩ解析における大気遅延によって生じる見かけの変位量に相当する。さらに実際のＧＰＳ解析では、２地点の観測局に置かれたＧＰＳ受信機内の時計誤差を相殺するために、異なる衛星についての位相差どうしをさらに引き算する。これを「二重位相差観測値」という。そこで、測位誤差シミュレーションを行う場合についても同様に、遅延量差の差を計算して見かけの変位量ベクトルＤを計算する。２つの異なる衛星Ｓ_１ とＳ_２ について遅延量の差の差（二重遅延量差） から行列Ｄを求めるときの行列Ｙは次のように書ける。
【００４１】
【数１１】

【００４２】
【数１２】

【００４３】
であり、計画行列Ａは、
【００４４】
【数１３】

【００４５】
と表される。より現実に近い測位誤差シミュレーションを行う場合は、４つないしそれ以上の衛星を想定して最小二乗法を実施する。すなわち、ｎ個の衛星について計算する場合は、行列ＹとＡの列の数が、一重差の場合はｎ×ｍ列となり、二重差の場合は（ｎ−１）×ｍ列となる。
【００４６】
また、先の式（８）で、仰角θの他に方位Ｒに依存した遅延勾配を考慮したマッピング関数ｍｆ（θ， α） （遅延勾配モデルとも言う）を大気遅延推定に用いたときのシミュレーションももちろん可能である。このときは、偏微分係数となるｍｆ（θ， α） は、大気の非一様性を単純な平面勾配で近似して構築され、ＭａｃＭｉｌｌａｎ（１９９６）を例に挙げると以下の式で表される。
【００４７】
【数１４】

【００４８】
右辺第１項が従来の球対称大気モデルのマッピング関数、右辺第２項が大気の非一様性の効果を示す。この第２項を求める手順としては、まず平面近似された大気勾配の大きさと向きをベクトルＧ（勾配ベクトル） で表す。上式で、Ｇの南北成分がＧ_Ｎ 、東西成分がＧ_Ｅ である。次に電波の到来方位αの単位ベクトルを考慮して大気遅延の水平勾配の強度と方位を表す量を求める。言い換えれば、この操作ではＧと電波到来方向のベクトルとの内積をとることになる。この第２項で乗じられているｃｏｔ（θ） は、仰角が減少するにしたがって水平方向の伝搬経路が長くなる結果、水平勾配の強度が強調される効果を考慮した項、及びｍ_０ （θ） は以上のすべての効果を天頂遅延に投影するための項である。ｍ_０ （θ） は従来の球対称マッピング関数を用いる。ここで、本来は第２項のｃｏｔ（θ） に含まれるθの大きさは、大気の非一様性を考慮したときの仰角θ_０ を用いるべきである。しかしながら、双方の差（θ’−θ） は仰角５度でも０．０３度程度に過ぎないため、近似的にθが使用されている。
【００４９】
ここでまとめると、遅延勾配モデルを評価対象として測位誤差シミュレーションを行う場合、天頂遅延量ＺＤの他に勾配ベクトルＧの南北成分Ｇ_Ｎ ，東西成分Ｇ_Ｅ 、の２つが推定パラメータとして加わることになる。
【００５０】
こうしてシミュレートされた測位誤差は、次のようにして実際の測位解の評価に使用できる。つまり、実際のＧＰＳやＶＬＳＩ測位がなされた同時期、同地域の非静力学モデルから計算した測位誤差と、実際の測位結果のばらつきと比べることで、測位結果の評価が可能になる。例えば、突発的な地殻変動に見えるような変位がＧＰＳによって解析された場合に、測位誤差シミュレーションで同様の変位が計算されれば、実際の地殻変動ではないと見なされるし、逆であれば地殻変動の可能性が高いと判断できる。また、実際の測位解の振る舞いと、シミュレートされた測位誤差との比較を様々な気象条件下で行うことによって、現行の大気モデルの改良に必要な基礎データを取得することができる。
【００５１】
図１は本発明の概念の説明図である。図１において、Ａは気象学の数値シミュレーションデータによる大気遅延量（視線遅延量）、測位誤差の評価システムを表す。Ｂは既存の解析システムであって、ＧＰＳもしくはＶＬＳＩによる測位観測システムである。本発明のシステムＡは、気象学に用いる数値シミュレーションデータ５０１をもとに数値計算５０２により、視線遅延量および測位誤差シミュレーションを行う。５０３は数値シミュレーションデータをもとにレイトレーシングによる視線遅延量計算の処理を表す。５０４は数値シミュレーションデータに基づいて５０３の処理で求めた視線遅延量により正確な測位誤差をシミュレーションにより求め、その測位誤差と実際の観測データにより測位観測結果と比較することにより実際の測位結果を評価する測位誤差シミュレーションの処理を表す。
【００５２】
既存の解析システムＢは、図１４で説明したものと同様である。球対称モデル、あるいは遅延勾配モデル等のモデルに基づく大気遅延量推定の処理４２３は、それぞれのモデルのマッピング関数により大気遅延量を推定する処理を表す（球対称モデルおよび遅延勾配モデルによる大気遅延量推定は図１５参照）。
【００５３】
図２は、本発明のシミュレーションによる測位誤差算出方法の概要を示す図である。図２において、６１０は、本発明で導入可能なベクトルＹの種別を表す。６０１は観測方程式を表す。６０２は観測量のベクトルＹを求める処理を表し、視線方向の遅延量の計算結果ΔＬ等である。６０３は未知パラメータのベクトルＤを表し、視線遅延量によって生じる見かけの変位（測位誤差）ベクトル、天頂遅延量ＺＤ等を要素とするものである。
【００５４】
６０４は計画行列Ａを表し、各パラメータの偏微分係数が要素であるものである。見かけの変位（シミュレーションにおける測位誤差）ベクトルは各成分の方向余弦である。天頂遅延量はＺＤであり、勾配遅延パラメータＧ_Ｎ 、Ｇ_Ｅ についてはマッピング関数で表す。６１２は最小二乗法によりベクトルＤを推定する処理を表す。
【００５５】
上記の観測方程式（１０）は既存の解析システムＢにおいても同じである。しかし、従来の測位誤差シミュレーションと本発明の測位誤差シミュレーションではそこで扱うパラメータが以下のように、行列Ｙ（ベクトルＹ）、行列Ｄ（ベクトルＤ）において相違する。
【００５６】
実際の測位において既存の解析システムＢにおける行列Ｙ、行列Ｄは次のようなものである。行列Ｙにおいて、エポック１番目からｍ番目まで、それぞれのエポックでのＧＰＳ衛星の位置（軌道情報Ｘ_Ｓ 、Ｙ_Ｓ 、Ｚ_Ｓ ）から地上の観測点（Ｘ_Ｒ 、Ｙ _Ｒ、Ｚ_Ｒ ） に到達した電波の伝搬時間に光速を乗じて計算した距離である。この量には、地上の観測点と衛星との間の真の距離情報、大気遅延量、衛星の時計誤差、観測点に設置した受信機の時計誤差などが含まれるため、真の距離と区別して疑似距離と呼ぶ。行列Ａは各推定要素の偏微分係数からなる計画行列である。
【００５７】
行列Ｄの未知パラメータにおいて、観測点位置（ｘ、ｙ、ｚ） 、及び大気による天頂方向の天頂遅延量（ＺＤ）は次のようなものである。既存の解析システムＢでは、その推定されるＺＤはあくまで、最小二乗法によってつじつまが合うように数学的に求められる値であり、実際の大気の状態がどのようになっているかとは無関係なものである。また、既存の解析システムＢでは大気中の水蒸気の変動が激しいときなどは、このＺＤの推定精度は落ちるものである。
【００５８】
本発明における測位誤差シミュレーションにおける行列Ｙ、行列Ｄは次のようなものである。行列Ｙは、エポック１番目からｍ番目まで、それぞれのエポックでのＧＰＳ衛星の位置（方位・仰角）をソフトウェアに入力し、数値シミュレーションデータを用いてレイトレーシングで計算した視線方向の大気遅延量である。したがって、本発明の測位誤差シミュレーションでは数値シミュレーションデータから求めた正確な大気状態に基づく伝搬遅延量が計算される。
【００５９】
本発明の測位誤差シミュレーションにおける行列Ｄの測位誤差ベクトル、及び大気による天頂方向の天頂遅延量（ＺＤ）は以下のようなものである。測位誤差ベクトルは、未知パラメータである観測点の真の位置（Ｘ_Ｒ 、Ｙ_Ｒ 、Ｚ_Ｒ ）から大気遅延によってずれた位置（Ｘ_Ａ 、Ｙ_Ａ 、Ｚ_Ａ ）までのベクトルである。すなわち、本発明の測位誤差シミュレーションにおいては、数値シミュレーションデータから求めた大気のみによって生じる伝搬遅延量が計算され、入力されているので、ここでの測位誤差ベクトルや天頂遅延量ＺＤは、明らかに大気によって生じた量と見なすことができる。したがって、本発明の測位誤差シミュレーションでは、推定される測位誤差や天頂遅延量は真の値に近いものである。また、本発明で計算した測位を使用すれば、入力する数値シミュレーションデータそのものの精度も向上し、より現実の大気の状態を推定することが可能になる。
【００６０】
さらに、実際の測位観測結果と本発明の測位誤差シミュレーションで推定された測位誤差とを比較して、双方の変動のパターンがほぼ同じであれば、実際の測位観測結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測結果の変動が地殻変動によるものと判定できる。
【００６１】
上記のように、本発明によれば、気象学の数値シミュレーションデータに基づいて正確に視線遅延量を計算し、その視線遅延量に基づいて測位誤差シミュレーションを行うことができる。そのため正確な測位を容易に把握することができるようになる。また、実際の観測で得た測位観測結果との相関や、そのデータを時系列データとして配列してその傾向を見ることが容易になり、測位観測結果のバラツキ等の原因を容易に判断できるようになる。また、実際の観測データをもとに測位計算する解析システムで採用したマッピング関数の精度、測位誤差の評価等も容易に行うことができる。
【００６２】
【発明の実施の形態】
図３は本発明の測位誤差シミュレーションシステムの構成図である。図３において、１はＣＰＵであり、測位誤差シミュレーションの演算を行うものである。２はメモリであり、測位誤差演算において必要なデータを記憶するものである。３は入出力インターフェイスである。４はデータ入力装置である、本発明で使用する気象学の数値シミュレーションデータ等を入力するものであり、６はデータベース作成手段であって、本発明で使用する数値シミュレーションデータ、本発明の測位誤差シミュレーションで求めた測位誤差等をデータベース化するものである。７は測位誤差シミュレーション手段であり、視線遅延量計算手段１１、測位誤差計算手段１２、評価手段１３を備えるものである。８はデータ保存装置であり、数値シミュレーションデータデータベース２１、マッピング関数保持手段２２、測位誤差シミュレーション結果データベース２３、測位観測結果データベース２４、各種データ保持手段２５を備えるものである。
【００６３】
図４は本発明の測位誤差シミュレーションシステムの装置構成の詳細である。図４の構成において、視線遅延量計算手段１１は前述の式（１）、（２）、（３）と（４）の数式の演算手段を備え、レイトレーシング法によりＧＰＳ方向の視線遅延量および天頂遅延量を求めるものである。測位誤差計算手段１２は前述の行列Ｙ、行列Ａ、行列Ｄにより最小二乗法により観測方程式の演算を行うものである。測位観測システム５１は既存のもので、図１の解析システムＢと同じである。
【００６４】
図４の構成において、数値シミュレーションデータ選択手段３４は、他システム３１から数値シミュレーションデータを入力する。数値シミュレーションデータは、例えば、非静力学モデルで計算されたものであり、大気の状態を正確に表しているもので、例えば、インターネット等のネットワークを介して入力する。あるいは、記憶媒体に格納されている状態で入力することもできる。もしくは、本発明のシステムに備える数値シミュレーションデータデータベース２１から入力する。入力された数値シミュレーションデータはメモリ領域３７に保持される。視線遅延量計算手段１１は数値シミュレーションデータとレイトレーシング法によりＧＰＳ方向の視線遅延量および天頂遅延量を計算する。計算された視線遅延量は視線遅延量保持手段４１に保持される。測位誤差計算手段１２は、数値シミュレーションデータをもとに求めた視線遅延量をもとに、行列Ｙ、行列Ａ、行列Ｄにより観測方程式の演算を行い、測位誤差を求める。そのようにして求めた測位誤差は測位誤差シミュレーション結果保持手段４２に保持される。評価手段１３は測位誤差シミュレーション結果の測位誤差と、既存の測位観測システム５１で得た測位誤差とを入力し、比較する。その統計的処理を行い、既存の測位観測誤差のバラツキ等の原因を視覚的に評価できるようにする。出力手段５４は評価手段１３で求めた評価結果をディスプレイに表示する、もしくはプリンタに印刷出力する。あるいは、記憶媒体に保存する。
【００６５】
入力した数値シミュレーションデータはデータベース作成手段３８により、数値シミュレーションデータ情報とともに数値シミュレーションデータデータベース２１および数値シミュレーションデータ情報保持手段３５に保存する。本発明の測位誤差シミュレーションでもとめた測位誤差は測位誤差シミュレーション結果データベース２３に保存する。また、既存の測位観測システム５１で得た測位誤差（測位観測データ）は測位観測結果データベース２４に保存することもできる。
【００６６】
図５は本発明の視線遅延量算出のための動作説明図である。本発明の視線遅延量算出は、数値シミュレーションデータとＧＰＳ衛星の位置に基づいて前述の式（１）、（２）、（３）の計算を行うものである。レイトレーシング法により視線方向（ＧＰＳ方向）の遅延量を求める。また天頂遅延量をもとめる。天頂遅延量は、数値シミュレーションデータをもとに計算する方法と視線遅延量をもとにマッピング関数を使用して求める等さまざまな方法で計算できる。
【００６７】
図５の構成の動作を説明する。視線遅延量計算手段１１は、数値シミュレーションデータ入力処理１１１により数値シミュレーションデータを入力する。数値シミュレーションデータとしては、気温、気圧、水蒸気量（露点温度）であり、各データは位置を示す格子点、年月日、時刻をもつ。入力された数値シミュレーションデータはメモリ領域１２０に保持される。屈折率算出手段１１２において、静水圧項計算手段１１３により静水圧項を計算する。湿潤項計算手段１１４は数値シミュレーションデータの気温、気圧、水蒸気量（露点温度）に基づいて、式（２）により湿潤項を計算する。静水圧項遅延量算出手段１１５と湿潤遅延量算出手段１１６により式（１）の第一項の計算をし、静水圧遅延量と湿潤遅延量をメモリ領域１２０に保持する。また、（Ｓ−Ｇ）計算手段１１７は数値シミュレーションデータおよびＧＰＳ衛星位置から式（１）の（Ｓ−Ｇ）の計算をし、メモリ領域１２０に保持する。静水圧遅延量と湿潤遅延量および（Ｓ−Ｇ）により視線方向遅延量ΔＬを求める。またレイトレーシング法により天頂遅延量を求める。天頂遅延量は視線遅延量をもとにマッピング関数で求める方法もある。使用するマッピング関数はマッピング関数指定入力１３１により指定され、選択される。算出した視線遅延量は視線遅延量保持手段４１に保持する（天頂遅延量の保持データは図示を省略している）。
【００６８】
図６は本発明の測位誤差算出のための動作説明図である。本発明の測位誤差シミュレーションは数値シミュレーションデータにより計算した視線遅延量をもとに、ＶＬＳＩでは前述の式（５）、（６）、（７）、（８）、（９）、（１０）により測位誤差シミュレーションを行う。ＧＰＳでは、前述の式（５）、（６）、（９）、（１０）、（１１）、（１２）、（１３）により測位誤差シミュレーションを行う。
【００６９】
一重位相差により測位誤差シミュレーションを行うか、二重位相差による測位誤差シミュレーションを行うか決め、また使用する大気モデルを決めて、測位誤差計算手段１２に数値シミュレーションデータで求めた各エポックの衛星方向の視線遅延量を入力し、行列Ｙを計算する（ベクトルＹ演算１５１の処理）。求めた行列Ｙはメモリ領域１６０に保持する。また、計画行列Ａを生成し（行列Ａ計算１５３の処理）、メモリ領域１６０に保持する。行列Ｄ（位置座標ｘ、ｙ、ｚ、天頂遅延量ＺＤが未知パラメータ）を設定し（未知パラメータＤ設定の処理１５２）、メモリ領域１６０に保持する。最小二乗法により誤差ベクトルｅが最小になるように、観測方程式（１０）に従って行列Ｄ（ ｘ、ｙ、ｚ、ＺＤ） を求める。このようにして求めたｘ、ｙ、ｚは測位誤差であり、ＺＤは推定天頂遅延量である。求めた測位誤差、推定天頂遅延量は測位誤差シミュレーション結果保持手段１８０に保持する。求めた測位誤差シミュレーションで求めた測位誤差はデータベース作成手段３８により測位誤差シミュレーション結果データベース２３に保存される。
【００７０】
測位誤差を評価する場合には、本発明に測位誤差シミュレーションで求めた測位誤差と実際の測位観測結果を評価手段１３に入力する。評価手段１３は測位誤差シミュレーションにより推定された測位誤差とを比較し、相違を視覚的にわかるように出力する。統計的に処理した双方の変動のパターンがほぼ同じであれば、実際の測位解析結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測システムで得た測位観測結果の変動が地殻変動によるものであると判定できる。
【００７１】
図６において、実際の測位観測システムから測位観測結果を入力する場合には、測位観測データ情報保持手段１７１を参照し、選択された測位観測データ情報選択手段１７２により測位観測データを選択して、測位観測結果データベースから入力する。
【００７２】
図７は、既存の解析システムで観測した測位観測結果をデータベースにするための構成である。測位観測結果入力手段２０１により測位観測システム２００から測位観測結果を入力する。データベース作成手段３８により、測位観測結果データベース情報を作成し、測位観測結果情報保持手段５２に保持するとともに、測位観測結果は測位観測結果データベース２４に保存する。
【００７３】
図８は本発明の視線遅延量算出のフローチャートである。図４、図５を参照し、図８の説明をする。数値シミュレーションデータ選択手段３４で選択された数値シミュレーションデータが入力され、数値シミュレーションデータを作成したモデル、数値シミュレーションデータの対象とする年月日、時刻、位置等を選択する（Ｓ１）。選択した数値シミュレーションデータを入力し、保持する（Ｓ２、Ｓ３）。選択した数値シミュレーションデータをデータベースに登録する場合には、データベース登録処理をし、データベースに登録する（Ｓ４、Ｓ５、Ｓ６、Ｓ７）。Ｓ５でデータベース登録の処理を行わないのであれば、数値シミュレーションデータから気圧、気温、水蒸気分布を求める（Ｓ８）。静水圧項の係数Ｋ_１ 、湿潤項の係数Ｋ_２ を定める（Ｓ９）。静水圧項Ｋ_１ （Ｐ／Ｔ）と湿潤項Ｋ_２ （Ｐｕ／Ｔ^２ ）を計算する（Ｓ１０）。屈折率に関する式（ｎ−１）×１０^６ を計算する（Ｓ１１）。（Ｓ−Ｇ）の計算をする（Ｓ１１’）。静水圧遅延ΔＬ_ｈ と湿潤遅延ΔＬ_ｗ 、それぞれの天頂遅延量ΔＬ_ｈ ^ｚ 、ΔＬ_ｗ ^ｚ を計算する（Ｓ１２）。静水圧遅延ΔＬ_ｈ と湿潤遅延ΔＬ_ｗ 、天頂遅延量ΔＬ_ｈ ^ｚ 、ΔＬ_ｗ ^ｚ を保持する（Ｓ１３）。
【００７４】
図９は本発明の測位誤差算出のフローチャートである。図９は一重位相差の場合を示す。図９の説明で図６を参照する。まず、ΔＬ_ｈ とΔＬ_ｗ を入力する（Ｓ１）。各イベント（各位置、視線方向）で、式（５）、もしくは式（１２）を計算する。ベクトルＹ（行列Ｙ）を計算する（Ｓ３）。計画行列Ａを求める（Ｓ４）。ベクトルＤを設定する（Ｓ５）。最小二乗法により誤差ベクトルｅを最小にするベクトルＤのパラメータ値を計算する（Ｓ６）。計算したベクトルＤの行列要素を測位誤差として保存する（Ｓ７）。
【００７５】
図１０は本発明の評価手段の視線遅延量の評価のフローチャートである。評価手段は、本発明の測位誤差シミュレーションシステムで求めた視線遅延量と、実際の測位観測システムでマッピング関数をもとに推定した視線遅延量を入力する（Ｓ１）。実際の測位観測システムで推定し視線遅延量と測位誤差シミュレーションで算出した視線遅延量を比較する（Ｓ２）。計算された誤差の可視化、統計処理を行う（Ｓ３）。可視化されたデータや統計処理結果を保存する（Ｓ４）。
【００７６】
図１１は本発明の評価手段の測位誤差の評価のフローチャートである。測位誤差シミュレーションで求めた測位誤差を入力する（Ｓ１）。測位誤差シミュレーションで求めた結果と実際の測位観測システムの観測結果を比較し、双方の相関や時系列のデータを可視化する（Ｓ２）。可視化されたデータや時系列データを保存する（Ｓ３）。
【００７７】
図１２は、本発明の各種データの構成を示す。（ａ）は数値シミュレーションデータの例を示す。（ｂ）は本発明における測位誤差シミュレーションにおける視線遅延量のデータ構成の例を示す。各ＧＰＳ衛星についてもつものである。（ｃ）は本発明における測位誤差シミュレーションによる測位誤差の例を示す。各ＧＰＳ衛星についてもつものである。（ｄ）は実際の測位観測システムにおける測位観測データの例を示す。
【００７８】
【発明の効果】
本発明の測位誤差シミュレーションにおいては、数値シミュレーションデータから求めた大気のみによって生じる伝搬遅延量を計算することができる。したがって、本発明の測位誤差シミュレーションでは、大気のみにより推定される測位誤差や大気に依存しない天頂遅延量を得ることができる。また、本発明では計算した視線遅延量を用いることにより、高精度のマッピング関数を考案できる。
【００７９】
さらに、実際の測位観測システムにより得られた測位観測結果と本発明の測位誤差シミュレーションで推定された測位誤差とを比較したとき、双方の変動のパターンがほぼ同じであれば、実際の測位観測結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測結果の変動が地殻変動によるものと判定できる。
【００８０】
このように、本発明によれば、正確な測位を容易に把握することができるようになる。また、実際の観測で得た測位観測結果との相関や、そのデータを時系列データとして配列してその傾向を見ることにより、測位観測結果のバラツキ等の原因を容易に判断できるようになる。また、実際の観測データをもとに測位計算する解析システムで採用したマッピング関数の精度、測位誤差評価等も容易に行うことができる。
【図面の簡単な説明】
【図１】本発明の概念の説明図である。
【図２】本発明の測位誤差算出方法の概要を示す図である。
【図３】本発明の測位誤差シミュレーションシステムの構成図である。
【図４】本発明の測位誤差シミュレーションシステムの装置構成の詳細を示す図である。
【図５】本発明の視線遅延量算出のための動作説明図である。
【図６】本発明の測位誤差算出のための動作説明図を示す。
【図７】既存の測位観測システムの測位観測結果についての本発明のデータベース作成の動作説明図を示す。
【図８】本発明の視線遅延量算出のフローチャートを示す図である。
【図９】本発明の測位誤差算出のフローチャートを示す図である。
【図１０】本発明の視線遅延量の評価のためのフローチャートを示す図である。
【図１１】本発明の測位誤差評価手段のフローチャートを示す図である。
【図１２】本発明の各種データの構成を示す図である。
【図１３】電波の大気伝搬遅延とマッピング関数の説明図を示す。
【図１４】既存の測位観測システムを示す図である。
【図１５】大気モデルの模式図を示す。
【符号の説明】
Ａ：大気遅延量・測位誤差評価システム
Ｂ：既存の解析システム
４０１：観測データ
４０２：測位計算システム
４１０：測位計算
４２３：大気遅延量推定
４２６：測位結果
５０１；気象学に用いる数値シミュレーションデータ
５０２：数値計算
５０３：レイトレーシングによる視線遅延量計算
５０４：実際のＧＰＳやＶＬＳＩ解析を考慮した測位誤差シミュレーション[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a positioning error simulation system, a positioning error evaluation method, a program for performing a positioning error simulation, and a storage medium for storing the program, and particularly to a positioning error simulation using VLSI or GPS. Calculate an accurate line-of-sight delay amount in the GPS satellite direction based on numerical simulation data that almost accurately represents the state of the atmosphere, and perform a positioning error simulation based on the line-of-sight delay amount to obtain an accurate positioning error. Can be. The numerical simulation data is obtained by numerical simulation of the state of the atmosphere as numerical data. In addition, based on the positioning error obtained by the positioning error simulation, the positioning error of the actual positioning observation result can be easily evaluated by comparing it with the positioning observation result obtained based on the actual observation data. .
[0002]
[Prior art]
Conventionally, position measurement uses radio waves transmitted from artificial satellites or quasars, and in the case of a VLSI, the position is measured based on the difference in the observation time of radio waves at two bases. In the case of the GPS, the distance to the satellite is estimated from the arrival time of the radio wave transmitted from the artificial satellite for a plurality of artificial satellites, and the position is measured.
[0003]
Radio waves transmitted from artificial satellites or quasars have a reduced propagation speed due to the refractive index in the atmosphere, and the propagation path has a curvature, which is apparently longer than the propagation path in a vacuum, and causes propagation delay in the atmosphere. FIG. 13A is an explanatory diagram of a propagation time delay (hereinafter, referred to as a line-of-sight delay amount) that occurs when a radio wave radiated from the GPS propagates in the atmosphere when positioning is performed by the GPS. G is a linear distance from a radio wave source such as a GPS satellite or a quasar to an antenna on the ground. S represents the actual path caused by atmospheric refraction. FIG. 13B is an explanatory diagram of a method of estimating a line-of-sight delay amount based on a zenith delay amount and a mapping function (described later).
[0004]
FIG. 14 shows a conventional positioning observation system. In FIG. 14, reference numeral 401 denotes observation data, such as a radio wave propagation time from a GPS satellite to a measurement point P, a GPS satellite azimuth, an elevation angle, and the like. Reference numeral 402 denotes a positioning calculation system which calculates a positioning error by positioning observation based on observation data. 410 is a positioning calculation means. Reference numeral 412 denotes a process of estimating a zenith delay amount, which is a process of estimating a zenith delay amount using observation data (radio wave propagation time) and a mapping function. Reference numeral 413 denotes a positioning error calculation process for calculating a positioning error based on the line-of-sight delay amount (zenith delay amount). 420 is a mapping function. Reference numeral 421 denotes the output of the positioning observation result.
[0005]
Conventionally, positioning error observations performed using observation data are performed by calculating a zenith delay amount based on the observed propagation time and calculating a positioning error based on the line-of-sight delay amount. Was.
[0006]
[Problems to be solved by the invention]
As described above, conventional positioning error observation is performed based on actual observation data that propagates in the atmosphere, and because of the propagation speed delay due to the refractive index in atmospheric propagation, the error due to the curvature generated in the propagation path, etc. Is observed differently each time measurement is performed due to the influence of the atmosphere. It was difficult to determine whether the change was due to atmospheric conditions or whether the observation position was actually moved due to a core change or the like.
[0007]
The present invention provides a positioning error simulation system that can calculate an accurate gaze delay based on numerical simulation data that accurately represents the state of the atmosphere, and that can obtain an accurate positioning error. It is an object of the present invention to be able to easily evaluate a positioning error observed in actual observation data by comparing with the above.
[0008]
[Means for Solving the Problems]
Numerical simulation data in meteorology divides the atmospheric space into a three-dimensional cubic lattice space, and calculates the state of the atmosphere in each lattice based on weather observation data. Smaller grids (hundreds of meters to several kilometers) and analysis with non-hydrostatic models give a fairly accurate representation of the state of the atmosphere. The present invention accurately calculates a line-of-sight delay amount of a radio wave from a GPS satellite to an observation point based on such numerical meteorological simulation data, and performs a positioning error simulation based on the line-of-sight delay amount. Enabled to calculate the positioning error.
[0009]
Prior to describing the embodiments of the present invention, the basic principle of atmospheric refraction will be described. Atmospheric delay due to atmospheric refraction is caused by the following two physical effects. One is an effect that the microwaves are decelerated more than in a vacuum because they pass through the neutral atmosphere, which is a dielectric medium. Secondly, since the propagation path has a curvature, the path actually extends from a straight line. At this time, since the density of the atmosphere is higher in the lower layer, the propagation path is upwardly convex. The delay ΔL due to these effects is called EPD (Access Path Delay) and is expressed by the following equation.
[0010]
(Equation 1)

[0011]
Here, G is a linear distance from a radio wave source such as a GPS satellite or a quasar to an antenna on the ground, and n (s) is a refractive index at a point s on a propagation path L where a curvature is caused by atmospheric refraction. At this time, the length of S is obtained by integrating the minute portion ds on L. The first term on the right side is the effect of deceleration, and the actual radio wave propagates on L at a propagation speed slower by the refractive index of the neutral atmosphere. However, in the analysis of VLSI or GPS, it is considered that the light beam propagates on the same route at a speed in a vacuum. Therefore, the propagation route is apparently extended by an amount obtained by multiplying the difference between the two speeds by the propagation time. The remaining second term on the right side has the effect of the curvature of the propagation path. Naturally, the delay due to this effect is zero in the zenith direction.
[0012]
Here, the calculation of the refractive index of the humid atmosphere can be obtained by the following equation (some steps in the middle are omitted here).
[0013]
(Equation 2)

[0014]
Here, P in the equation (2) is the total atmospheric pressure (hPa), Pv and T are the steam partial pressure (hPw) and the absolute temperature (K), respectively. Coefficient K of each term₁, K₂It is a constant determined from laboratory experiments and varies slightly between researchers. Generally, the first term on the right side of the equation (2) is referred to as “hydrostatic pressure term”. This term is a contribution proportional to the density of atmospheric components in hydrostatic equilibrium, and integrating this from the ground to the upper atmosphere gives a hydrostatic pressure delay. The second term on the right side of equation (2) is called a "wet term", and a delay obtained by integrating the term is called a wet delay.
[0015]
In the ray tracing calculation using the non-static model, the refractive index of each grid point is calculated using the equation (2) from the temperature, atmospheric pressure, and water vapor partial pressure given to each grid point of the non-static model. calculate. Next, the refraction angle is sequentially calculated from the calculated refractive index according to Snell's law on the path of the radio wave emitted from a certain observation point on the ground at an arbitrary elevation angle and direction to the satellite. This calculation itself is ray tracing. The conventional model can be evaluated by comparing the amount of delay in the gaze direction calculated by ray tracing with a conventional atmospheric model (this is called a mapping function). Next, a positioning error simulation according to the present invention will be described. First, a mapping function, which is a conventional atmospheric estimation method, will be briefly described, and then a positioning error simulation based on a result of ray tracing will be described.
[0016]
In the conventional atmospheric delay estimation, that is, in the analysis of the mapping function GPS or VLSI, the relative position between the observation stations included in the phase difference data of the radio wave between the two observation stations, the ionospheric delay, the atmospheric delay, the observation station or the artificial satellite An enormous number of observation equations are established using contributions such as clock errors in the system as unknown parameters, and each parameter is estimated by the least squares method. Speaking of atmospheric refraction, strictly speaking, there are temporal and spatial fluctuations in temperature, atmospheric pressure, and water vapor partial pressure, so that the atmospheric delay in the line of sight direction is a function of both time and space. However, if these are all unknowns, it is impossible to solve the observation equation because there are too many parameters. Therefore, it is assumed that the atmosphere maintains a uniform structure in the horizontal direction within a certain time, and during that time, the unknowns can be significantly reduced by performing modeling in which the atmospheric delay depends only on the elevation angle. Modeling also makes it easier to handle atmospheric delays when numerically calculating digitized observation data. A model called a mapping function was introduced based on this concept. Using the mapping function, the total atmospheric delay ΔL for a given elevation angle θ can be simply expressed as:
[0017]
(Equation 3)

[0018]
Where ΔL^z _hAnd ΔL^z _wAre the hydrostatic pressure delay and the wetting delay in the zenith direction, respectively. M multiplied by these_h(Θ) and M_w(Θ) are mapping functions for hydrostatic pressure delay and wet delay, respectively.
[0019]
As will be mentioned later, the mapping function is the partial derivative of atmospheric delay in the design matrix created when the observation equation is linearized. However, in VLSI, the mapping function difference between two points (“atmospheric delay difference” included in “single phase difference observation value” described later), and in GPS, the mapping function between two points obtained for different satellites The difference is a value obtained by taking the difference again between the differences (“difference in atmospheric delay difference” included in the “dual phase difference observed value” described later).
[0020]
The mapping function is basically the reciprocal of the sine of the elevation angle θ, that is, (1 / sin (θ)). However, when observing at a low elevation angle of 10 degrees or less, the curvature of the propagation path becomes tight because the microwave follows a longer path in the atmosphere. As a result, the delay due to the curvature effect increases exponentially. Moreover, the influence of the curvature of the earth cannot be ignored. Therefore, instead of a simple "reciprocal of the sine of the elevation angle", a mapping function in which the denominator is extended by a continued fraction as shown in the following equation (4) is considered as a basic form.
[0021]
(Equation 4)

[0022]
Here, the coefficients (a, b, c) of each term are constants. These coefficients are determined by comparing the atmospheric delay obtained by using ray tracing from the refractive index profile obtained from the spherically symmetric atmospheric model with the mapping function, and using the least-squares method so that both coincide. In the development of the actual mapping function, the coefficients are constants, variables that change with the atmospheric pressure and temperature of the ground weather data, and the final b / (sin θ + c) part is replaced with other function types such as tan θ Various improvements have been made by researchers to find more accurate mapping functions. The magnitude of the difference between the atmospheric delay amounts calculated from the mapping function developed in this way and the ray tracing method indicates the applicable range of the mapping function, that is, the estimation accuracy.
[0023]
As a specific procedure, a positioning error simulation caused by an atmospheric delay can be performed by inputting a ray-traced delay amount into an observation equation used in GPS or VLSI analysis and solving the observation equation.
[0024]
Specifically, it is calculated as follows. First, in the calculation of GPS or VLSI, a phase difference between two points is calculated, so that a difference between observation points is also calculated for a delay amount. Observation point R₁And R₂Assuming that radio waves arrive in parallel from a common satellite, the delay difference between the i-th two points in the observation epoch is
[0025]
(Equation 5)

[0026]
Where ΔL^S(R₁)_iAnd ΔL^S(R₂)_iIs the observation point R₁And R₂Is the amount of line-of-sight delay at an arbitrary azimuth and elevation. The suffixes S and i indicate the satellite and epoch numbers, respectively (the epoch is a delimiter for observation). Here, the observation point R₁Is a fixed point, the delay amount difference ΔL^S(R₂-R₁)_iIs a delay amount difference between two points, and corresponds to an atmospheric delay amount difference included in a “single phase difference observed value” in a normal GPS analysis.
[0027]
Here, using the delay amount difference for each of a large number of epochs as observation data, a vector most suitable for this delay amount difference can be calculated by the least square method. That is, when this relationship is linearized,
[0028]
(Equation 6)

[0029]
It is. In the equation (6), the matrix e is an error vector, and the matrix Y is a vector of observation amounts from the first to the m-th epoch obtained from the delay residual.
[0030]
(Equation 7)

[0031]
It is expressed as The matrix A is a planning matrix composed of partial differential coefficients of each element. For simplicity, parameters to be estimated include an error vector of an observation point position (apparent “displacement vector” described later) and an atmospheric delay amount in a zenith direction. Then
[0032]
(Equation 8)

[0033]
Becomes here,
α_i, Β_iAnd γ_i(I = 1, 2,... M) and the like are the direction cosine of each of the x, y, and z components of the delay amount difference. Also, Δmf (θ) (R₂-R₁) Is the partial differential coefficient of the delay amount (here, the difference between two points of the mapping function with θ as the elevation angle). Further, D is a vector obtained as a result of fitting this equation,
[0034]
(Equation 9)

[0035]
Can be expressed as Here, x, y, and z are components of the vector, and are assumed to be north in the positive direction of the x-axis, east in the positive direction of the y-axis, and vertically upward in the positive direction of the z-axis. ZD is the amount of delay in the zenith direction.
[0036]
This vector D is calculated at the observation point R₂Is a vector that shifts the position from the true position, and can be regarded as an apparent displacement vector, that is, a positioning error caused by atmospheric delay. The calculation of this apparent displacement vector is a positioning error simulation, and the behavior of the positioning error can be examined by examining the magnitude and direction of this vector under various atmospheric conditions.
[0037]
Considering that the apparent displacement vector D is actually obtained from the delay difference using the least squares method, the observation equation in this case is as follows.
[0038]
(Equation 10)

[0039]
Can be written as
[0040]
The matrix D obtained here corresponds to, for example, an apparent displacement amount caused by atmospheric delay in VLSI analysis. Furthermore, in the actual GPS analysis, the phase difference between different satellites is further subtracted in order to cancel the clock error in the GPS receivers located at the two observation stations. This is called “double phase difference observation value”. Therefore, also in the case where the positioning error simulation is performed, the difference of the delay amount difference is calculated to calculate the apparent displacement vector D. Two different satellites S₁And S₂The matrix Y for obtaining the matrix D from the difference between the delay amounts (double delay amount difference) can be written as follows.
[0041]
[Equation 11]

[0042]
(Equation 12)

[0043]
And the design matrix A is
[0044]
(Equation 13)

[0045]
It is expressed as When a more realistic positioning error simulation is performed, the least squares method is performed assuming four or more satellites. That is, when calculating for n satellites, the number of columns of the matrices Y and A is n × m columns for a single difference, and (n−1) × m columns for a double difference.
[0046]
In addition, in the above equation (8), a simulation is performed when a mapping function mf (θ, α) (also referred to as a delay gradient model) that takes into account a delay gradient dependent on the azimuth R in addition to the elevation angle θ is used for atmospheric delay estimation. Of course, it is possible. In this case, mf (θ, α), which is a partial differential coefficient, is constructed by approximating the non-uniformity of the atmosphere with a simple plane gradient, and is expressed by the following equation using MacMillan (1996) as an example. You.
[0047]
[Equation 14]

[0048]
The first term on the right-hand side shows the mapping function of the conventional spherically symmetric atmosphere model, and the second term on the right-hand side shows the effect of atmospheric non-uniformity. As a procedure for obtaining the second term, first, the magnitude and direction of the plane-graded atmospheric gradient are represented by a vector G (gradient vector). In the above equation, the north-south component of G is G_NAnd the east-west component is G_EIt is. Next, an amount representing the strength and direction of the horizontal gradient of the atmospheric delay is determined in consideration of the unit vector of the arrival direction α of the radio wave. In other words, in this operation, the inner product of G and the vector of the radio wave arrival direction is obtained. The term cot (θ) multiplied by the second term is a term that takes into account the effect that the intensity of the horizontal gradient is enhanced as a result that the horizontal propagation path becomes longer as the elevation angle decreases, and m₀(Θ) is a term for projecting all the above effects to the zenith delay. m₀(Θ) uses a conventional spherically symmetric mapping function. Here, the magnitude of θ originally contained in the second term cot (θ) is the elevation angle θ in consideration of the non-uniformity of the atmosphere.₀Should be used. However, the difference ([theta] '-[theta]) between them is only about 0.03 [deg.] Even when the elevation angle is 5 [deg.], So [theta] is approximately used.
[0049]
In summary, when the positioning error simulation is performed with the delay gradient model as the evaluation target, the north-south component G of the gradient vector G in addition to the zenith delay amount ZD_N, East-West component G_E, Are added as estimation parameters.
[0050]
The positioning error simulated in this way can be used for evaluating an actual positioning solution as follows. That is, at the same time when the actual GPS or VLSI positioning is performed, the positioning result calculated by comparing the positioning error calculated from the non-static model in the same area with the dispersion of the actual positioning result can be evaluated. For example, when a displacement that appears to be sudden crustal deformation is analyzed by GPS, if a similar displacement is calculated by positioning error simulation, it is regarded as not an actual crustal deformation, and if it is the opposite, the crustal deformation is considered. It can be determined that the possibility of fluctuation is high. In addition, by comparing the behavior of the actual positioning solution with the simulated positioning error under various weather conditions, it is possible to obtain basic data necessary for improving the current atmospheric model.
[0051]
FIG. 1 is an explanatory diagram of the concept of the present invention. In FIG. 1, A represents an atmospheric delay amount (gaze delay amount) and a positioning error evaluation system based on numerical simulation data of meteorology. B is an existing analysis system, which is a positioning observation system using GPS or VLSI. The system A of the present invention performs a gaze delay amount and a positioning error simulation by numerical calculation 502 based on numerical simulation data 501 used for meteorology. Reference numeral 503 denotes a process of calculating a visual line delay amount by ray tracing based on numerical simulation data. Reference numeral 504 denotes an accurate positioning error obtained by simulation based on the visual line delay amount obtained in the process of 503 based on the numerical simulation data, and evaluates the actual positioning result by comparing the positioning error with the positioning observation result based on the actual observation data. Represents a positioning error simulation process.
[0052]
The existing analysis system B is the same as that described with reference to FIG. Atmospheric delay amount estimation processing 423 based on a spherical symmetric model or a delay gradient model or the like represents a process of estimating the atmospheric delay amount by a mapping function of each model (atmospheric delay amount by a spherical symmetric model and a delay gradient model). The estimation is shown in FIG. 15).
[0053]
FIG. 2 is a diagram showing an outline of a positioning error calculation method by simulation according to the present invention. In FIG. 2, reference numeral 610 denotes a type of a vector Y that can be introduced in the present invention. 601 represents an observation equation. Reference numeral 602 denotes a process for obtaining a vector Y of the observation amount, such as a calculation result ΔL of the delay amount in the line of sight. Reference numeral 603 denotes a vector D of an unknown parameter, and includes an apparent displacement (positioning error) vector caused by a line-of-sight delay amount, a zenith delay amount ZD, and the like.
[0054]
Reference numeral 604 denotes a design matrix A, in which the partial differential coefficient of each parameter is an element. The apparent displacement (positioning error in the simulation) vector is the direction cosine of each component. The zenith delay amount is ZD, and the gradient delay parameter G_N, G_EIs represented by a mapping function. Reference numeral 612 denotes processing for estimating the vector D by the least square method.
[0055]
The above observation equation (10) is the same in the existing analysis system B. However, the parameters handled in the conventional positioning error simulation and the positioning error simulation of the present invention are different between the matrix Y (vector Y) and the matrix D (vector D) as follows.
[0056]
In the actual positioning, the matrices Y and D in the existing analysis system B are as follows. In the matrix Y, the positions of the GPS satellites at each epoch (orbit information X_S, Y_S, Z_S) To the ground observation point (X_R, Y_R, Z_R) Is the distance calculated by multiplying the propagation time of the radio wave arriving at by the speed of light. This amount includes the true distance information between the observation point on the ground and the satellite, the atmospheric delay, the clock error of the satellite, and the clock error of the receiver installed at the observation point. Separately, it is called pseudo distance. The matrix A is a design matrix composed of partial derivatives of each estimation element.
[0057]
Among the unknown parameters of the matrix D, the observation point position (x, y, z) and the zenith delay (ZD) in the zenith direction due to the atmosphere are as follows. In the existing analysis system B, the estimated ZD is a value that is mathematically obtained so as to be consistent by the least squares method, and is irrelevant to the actual state of the atmosphere. It is. Further, in the existing analysis system B, when the water vapor in the atmosphere fluctuates greatly, the estimation accuracy of the ZD decreases.
[0058]
The matrices Y and D in the positioning error simulation according to the present invention are as follows. The matrix Y is the amount of atmospheric delay in the line-of-sight direction calculated by ray tracing using numerical simulation data by inputting the positions (azimuth and elevation) of the GPS satellites at each epoch from the first epoch to the m-th epoch. is there. Therefore, in the positioning error simulation of the present invention, the propagation delay amount based on the accurate atmospheric condition obtained from the numerical simulation data is calculated.
[0059]
The positioning error vector of the matrix D and the zenith delay (ZD) in the zenith direction due to the atmosphere in the positioning error simulation of the present invention are as follows. The positioning error vector is obtained by calculating the true position (X_R, Y_R, Z_R) Due to atmospheric delay (X_A, Y_A, Z_A). That is, in the positioning error simulation of the present invention, since the propagation delay amount caused only by the atmosphere obtained from the numerical simulation data is calculated and input, the positioning error vector and the zenith delay amount ZD here are obviously Can be considered as the amount produced by Therefore, in the positioning error simulation of the present invention, the estimated positioning error and the zenith delay amount are close to the true values. Further, if the positioning calculated by the present invention is used, the accuracy of the input numerical simulation data itself is also improved, and it is possible to estimate the actual state of the atmosphere.
[0060]
Furthermore, comparing the actual positioning observation result with the positioning error estimated by the positioning error simulation of the present invention, if the patterns of both fluctuations are substantially the same, the error due to the atmosphere in the actual positioning observation result is remarkable. It means that it is included, and if the patterns are not similar, it can be determined that the fluctuation of the actual positioning observation result is due to crustal deformation.
[0061]
As described above, according to the present invention, the gaze delay amount can be accurately calculated based on meteorological numerical simulation data, and a positioning error simulation can be performed based on the gaze delay amount. Therefore, accurate positioning can be easily grasped. In addition, it is easy to see the correlation with the positioning observation result obtained by actual observation and to see the tendency by arranging the data as time-series data, and to easily determine the cause of the variation in the positioning observation result. become. In addition, the accuracy of the mapping function and the evaluation of the positioning error used in the analysis system that performs the positioning calculation based on the actual observation data can be easily performed.
[0062]
BEST MODE FOR CARRYING OUT THE INVENTION
FIG. 3 is a configuration diagram of the positioning error simulation system of the present invention. In FIG. 3, reference numeral 1 denotes a CPU for performing calculation of a positioning error simulation. Reference numeral 2 denotes a memory for storing data necessary for positioning error calculation. 3 is an input / output interface. Numeral 4 denotes a data input device for inputting numerical simulation data of meteorology used in the present invention. Numeral 6 denotes a database creating means, which is a numerical simulation data used in the present invention and a positioning error of the present invention. A database is used for positioning errors and the like obtained by simulation. Reference numeral 7 denotes a positioning error simulation means, which includes a line-of-sight delay amount calculation means 11, a positioning error calculation means 12, and an evaluation means 13. A data storage device 8 includes a numerical simulation data database 21, a mapping function holding unit 22, a positioning error simulation result database 23, a positioning observation result database 24, and various data holding units 25.
[0063]
FIG. 4 shows the details of the device configuration of the positioning error simulation system of the present invention. In the configuration shown in FIG. 4, the line-of-sight delay amount calculating means 11 includes the calculating means of the above-described equations (1), (2), (3) and (4). The zenith delay amount is obtained. The positioning error calculation means 12 performs the calculation of the observation equation by the least square method using the aforementioned matrix Y, matrix A, and matrix D. The positioning observation system 51 is an existing system, and is the same as the analysis system B of FIG.
[0064]
In the configuration shown in FIG. 4, the numerical simulation data selection means 34 inputs numerical simulation data from another system 31. The numerical simulation data is, for example, data calculated by a non-static model, accurately represents the state of the atmosphere, and is input via a network such as the Internet, for example. Alternatively, the input can be made in a state of being stored in a storage medium. Alternatively, the data is input from the numerical simulation data database 21 provided in the system of the present invention. The input numerical simulation data is held in the memory area 37. The line-of-sight delay amount calculation means 11 calculates the line-of-sight delay amount and the zenith delay amount in the GPS direction using numerical simulation data and the ray tracing method. The calculated line-of-sight delay amount is held in the line-of-sight delay amount holding unit 41. The positioning error calculation means 12 calculates an observation equation by using the matrix Y, the matrix A, and the matrix D based on the line-of-sight delay amount obtained based on the numerical simulation data, and obtains a positioning error. The positioning error thus obtained is held in the positioning error simulation result holding means 42. The evaluation means 13 inputs and compares the positioning error of the positioning error simulation result with the positioning error obtained by the existing positioning observation system 51. By performing the statistical processing, it is possible to visually evaluate a cause such as the variation of the existing positioning observation error. The output unit 54 displays the evaluation result obtained by the evaluation unit 13 on a display or prints out the result on a printer. Alternatively, it is stored in a storage medium.
[0065]
The inputted numerical simulation data is stored in the numerical simulation data database 21 and the numerical simulation data information holding means 35 together with the numerical simulation data information by the database creating means 38. The positioning error determined by the positioning error simulation of the present invention is stored in the positioning error simulation result database 23. Further, positioning errors (positioning observation data) obtained by the existing positioning observation system 51 can be stored in the positioning observation result database 24.
[0066]
FIG. 5 is an explanatory diagram of the operation for calculating the gaze delay amount according to the present invention. In the gaze delay amount calculation of the present invention, the above-described equations (1), (2), and (3) are calculated based on the numerical simulation data and the position of the GPS satellite. The delay amount in the line-of-sight direction (GPS direction) is obtained by the ray tracing method. Also, determine the amount of zenith delay. The zenith delay amount can be calculated by various methods such as a method of calculating based on numerical simulation data and a method of obtaining the zenith delay amount using a mapping function based on the line-of-sight delay amount.
[0067]
The operation of the configuration of FIG. 5 will be described. The line-of-sight delay amount calculation means 11 inputs numerical simulation data by the numerical simulation data input processing 111. Numerical simulation data includes temperature, atmospheric pressure, and water vapor amount (dew point temperature), and each data has a grid point indicating a position, a date, and a time. The input numerical simulation data is held in the memory area 120. In the refractive index calculating means 112, the hydrostatic pressure term is calculated by the hydrostatic pressure calculating means 113. The wetting term calculation means 114 calculates the wetting term by equation (2) based on the temperature, pressure, and water vapor amount (dew point temperature) of the numerical simulation data. The first term of the equation (1) is calculated by the hydrostatic pressure term delay calculating means 115 and the wet delay calculating means 116, and the hydrostatic pressure delay and the wet delay are stored in the memory area 120. Further, the (SG) calculation means 117 calculates (SG) of Expression (1) from the numerical simulation data and the GPS satellite position, and stores the calculation in the memory area 120. The gaze direction delay amount ΔL is obtained from the hydrostatic pressure delay amount, the wet delay amount, and (SG). Further, the zenith delay amount is obtained by a ray tracing method. The zenith delay amount may be obtained by a mapping function based on the line-of-sight delay amount. The mapping function to be used is designated by the mapping function designation input 131 and is selected. The calculated line-of-sight delay amount is stored in the line-of-sight delay amount holding means 41 (data for holding the zenith delay amount is not shown).
[0068]
FIG. 6 is an explanatory diagram of an operation for calculating a positioning error according to the present invention. In the positioning error simulation of the present invention, the VLSI is based on the line-of-sight delay amount calculated from the numerical simulation data according to the above equations (5), (6), (7), (8), (9), and (10). Perform a positioning error simulation. In the GPS, a positioning error simulation is performed by the above equations (5), (6), (9), (10), (11), (12), and (13).
[0069]
Determine whether to perform the positioning error simulation using the single phase difference or the positioning error simulation using the double phase difference, determine the atmospheric model to be used, and instruct the positioning error calculation means 12 to obtain the satellite direction of each epoch obtained from the numerical simulation data. Is input, and the matrix Y is calculated (the processing of the vector Y operation 151). The obtained matrix Y is stored in the memory area 160. Further, a planning matrix A is generated (processing of the matrix A calculation 153), and stored in the memory area 160. A matrix D (position coordinates x, y, z, and zenith delay amount ZD are unknown parameters) is set (process 152 for setting unknown parameters D) and stored in the memory area 160. A matrix D (x, y, z, ZD) is obtained according to the observation equation (10) so that the error vector e is minimized by the least square method. X, y, and z thus obtained are positioning errors, and ZD is an estimated zenith delay amount. The obtained positioning error and the estimated zenith delay amount are stored in the positioning error simulation result storage unit 180. The positioning error obtained by the obtained positioning error simulation is stored in the positioning error simulation result database 23 by the database creating means 38.
[0070]
When evaluating the positioning error, the positioning error obtained by the positioning error simulation according to the present invention and the actual positioning observation result are input to the evaluation means 13. The evaluation means 13 compares the positioning error estimated by the positioning error simulation with the positioning error and outputs the difference so as to be visually recognized. If the patterns of both statistically processed fluctuations are almost the same, it means that the actual positioning analysis results include significant errors due to the atmosphere.If the patterns are not similar, the results can be obtained with the actual positioning observation system. It can be determined that the fluctuation of the positioning observation result is due to crustal deformation.
[0071]
In FIG. 6, when the positioning observation result is input from the actual positioning observation system, the positioning observation data is selected by the selected positioning observation data information selecting unit 172 with reference to the positioning observation data information holding unit 171. Input from the positioning observation result database.
[0072]
FIG. 7 shows a configuration for converting a positioning observation result observed by an existing analysis system into a database. The positioning observation result is input from the positioning observation system 200 by the positioning observation result input means 201. The database creation unit 38 creates positioning observation result database information and holds it in the positioning observation result information holding unit 52, and saves the positioning observation result in the positioning observation result database 24.
[0073]
FIG. 8 is a flowchart of the gaze delay amount calculation according to the present invention. FIG. 8 will be described with reference to FIGS. The numerical simulation data selected by the numerical simulation data selection means 34 is input, and the model that created the numerical simulation data, and the date, time, position, and the like to be subjected to the numerical simulation data are selected (S1). The selected numerical simulation data is input and held (S2, S3). When registering the selected numerical simulation data in the database, a database registration process is performed and registered in the database (S4, S5, S6, S7). If the database registration processing is not performed in S5, the atmospheric pressure, temperature, and water vapor distribution are obtained from the numerical simulation data (S8). Hydrostatic pressure coefficient K₁, The coefficient of the wet term K₂Is determined (S9). Hydrostatic pressure term K₁(P / T) and wet term K₂(Pu / T²) Is calculated (S10). Formula (n-1) × 10 relating to refractive index⁶Is calculated (S11). (S-G) is calculated (S11 '). Hydrostatic pressure delay ΔL_hAnd wet delay ΔL_w, Each zenith delay amount ΔL_h ^z, ΔL_w ^zIs calculated (S12). Hydrostatic pressure delay ΔL_hAnd wet delay ΔL_w, Zenith delay amount ΔL_h ^z, ΔL_w ^zIs held (S13).
[0074]
FIG. 9 is a flowchart of the positioning error calculation according to the present invention. FIG. 9 shows the case of a single phase difference. FIG. 6 will be referred to in the description of FIG. First, ΔL_hAnd ΔL_wIs input (S1). Formula (5) or Formula (12) is calculated for each event (each position, line-of-sight direction). A vector Y (matrix Y) is calculated (S3). A design matrix A is obtained (S4). A vector D is set (S5). The parameter value of the vector D that minimizes the error vector e is calculated by the least square method (S6). The calculated matrix element of the vector D is stored as a positioning error (S7).
[0075]
FIG. 10 is a flowchart of the evaluation of the gaze delay amount by the evaluation means of the present invention. The evaluation means inputs the gaze delay amount obtained by the positioning error simulation system of the present invention and the gaze delay amount estimated based on the mapping function by the actual positioning observation system (S1). The gaze delay estimated by the actual positioning observation system and the gaze delay calculated by the positioning error simulation are compared (S2). Visualization and statistical processing of the calculated error are performed (S3). The visualized data and the result of the statistical processing are stored (S4).
[0076]
FIG. 11 is a flowchart of the evaluation of the positioning error by the evaluation means of the present invention. The positioning error obtained by the positioning error simulation is input (S1). The result obtained by the positioning error simulation is compared with the observation result of the actual positioning observation system, and both correlations and time series data are visualized (S2). The visualized data and time-series data are stored (S3).
[0077]
FIG. 12 shows the configuration of various data of the present invention. (A) shows an example of numerical simulation data. (B) shows an example of the data configuration of the line-of-sight delay amount in the positioning error simulation according to the present invention. It has about each GPS satellite. (C) shows an example of the positioning error by the positioning error simulation in the present invention. It has about each GPS satellite. (D) shows an example of positioning observation data in an actual positioning observation system.
[0078]
【The invention's effect】
In the positioning error simulation of the present invention, it is possible to calculate the propagation delay amount caused only by the atmosphere obtained from the numerical simulation data. Therefore, in the positioning error simulation of the present invention, it is possible to obtain a positioning error estimated only from the atmosphere and a zenith delay amount independent of the atmosphere. In the present invention, a highly accurate mapping function can be devised by using the calculated line-of-sight delay amount.
[0079]
Furthermore, when the positioning observation result obtained by the actual positioning observation system is compared with the positioning error estimated by the positioning error simulation of the present invention, if the two fluctuation patterns are substantially the same, the actual positioning observation result is obtained. Means that errors due to the atmosphere are remarkably included in the data, and if the patterns are not similar, it can be determined that the change in the actual positioning observation result is due to crustal deformation.
[0080]
As described above, according to the present invention, accurate positioning can be easily grasped. In addition, by observing the correlation with the positioning observation result obtained by actual observation and arranging the data as time-series data and observing the tendency, it is possible to easily determine the cause of the variation in the positioning observation result. In addition, the accuracy of the mapping function, the positioning error evaluation, and the like employed in the analysis system that performs the positioning calculation based on the actual observation data can be easily performed.
[Brief description of the drawings]
FIG. 1 is an explanatory diagram of the concept of the present invention.
FIG. 2 is a diagram showing an outline of a positioning error calculation method of the present invention.
FIG. 3 is a configuration diagram of a positioning error simulation system of the present invention.
FIG. 4 is a diagram showing details of a device configuration of a positioning error simulation system of the present invention.
FIG. 5 is an explanatory diagram of an operation for calculating a visual line delay amount according to the present invention.
FIG. 6 is a diagram illustrating an operation for calculating a positioning error according to the present invention.
FIG. 7 is an explanatory diagram illustrating an operation of creating a database according to the present invention with respect to positioning observation results of an existing positioning observation system.
FIG. 8 is a diagram showing a flowchart of gaze delay amount calculation according to the present invention.
FIG. 9 is a diagram showing a flowchart of a positioning error calculation of the present invention.
FIG. 10 is a diagram showing a flowchart for evaluating a gaze delay amount according to the present invention.
FIG. 11 is a diagram showing a flowchart of a positioning error evaluation unit of the present invention.
FIG. 12 is a diagram showing a configuration of various data of the present invention.
FIG. 13 is an explanatory diagram of an atmospheric propagation delay of a radio wave and a mapping function.
FIG. 14 is a diagram showing an existing positioning observation system.
FIG. 15 shows a schematic diagram of an atmospheric model.
[Explanation of symbols]
A: Atmospheric delay / positioning error evaluation system
B: Existing analysis system
401: Observation data
402: Positioning calculation system
410: Positioning calculation
423: Atmospheric delay estimation
426: Positioning result
501: Numerical simulation data used for meteorology
502: Numerical calculation
503: Gaze delay calculation by ray tracing
504: Positioning error simulation considering actual GPS and VLSI analysis

Claims (13)

Means for inputting numerical simulation data representing the state of the atmosphere, gaze delay calculating means for calculating a gaze delay based on the numerical simulation data, and a positioning error based on the gaze delay calculated by the gaze delay calculating means. A positioning error simulation system, comprising: a positioning error calculation means for obtaining a positioning error by simulation. The positioning error simulation system according to claim 1, wherein the numerical simulation data is obtained based on a non-hydrostatic numerical forecast model. An input means of the positioning observation result measured based on the actual observation data is compared with the positioning error obtained by the positioning error simulation and the positioning error simulation, and the positioning error of the positioning observation result obtained from the actual observation data is evaluated. The positioning error simulation system according to claim 1, further comprising an evaluation unit. 2. The positioning error simulation system according to claim 1, further comprising a database creation unit, and holding a numerical simulation data database of the inputted numerical simulation data and numerical simulation data information including conditions for creating the numerical simulation data. A positioning error database that holds a positioning error obtained by the positioning error simulation as a database together with positioning error generation information including a model of numerical simulation data used for its generation and a model represented by a mapping function. 2. The positioning error simulation system according to 1. 2. The positioning error simulation system according to claim 1, wherein the input actual positioning observation result is used as a positioning observation result database and held together with the positioning observation result information. Numerical simulation data representing the state of the atmosphere is input, stored in a memory, a gaze delay amount is calculated by an arithmetic unit based on the numerical simulation data, the gaze delay amount is held, and the gaze delay amount is calculated by the arithmetic unit. The positioning error is calculated by the positioning error simulation based on the positioning error, the positioning error is held, the actual positioning observation result is input, and the positioning observation result is included in the positioning observation result by comparing with the result calculated by the positioning error simulation. A positioning error evaluation method, wherein a positioning error is evaluated. The positioning error evaluation method according to claim 7, wherein the numerical simulation data is generated based on a non-hydrostatic numerical forecast model. The positioning error evaluation method according to claim 7, wherein the input numerical simulation data is used as a numerical simulation data database together with numerical simulation data information including conditions for creating numerical simulation data. The positioning error evaluation method according to claim 7, wherein the positioning error is held as a database together with positioning error generation information including a model of numerical simulation data used in the positioning error simulation and a model represented by a mapping function. The positioning error evaluation method according to claim 7, wherein the input actual positioning observation result is stored as a positioning observation data database together with the positioning observation result information. Numerical simulation data representing the state of the atmosphere is input, stored in a memory, a gaze delay amount is calculated by an arithmetic unit based on the numerical simulation data, the gaze delay amount is held, and the gaze delay amount is calculated by the arithmetic unit. A positioning error is calculated by a positioning error simulation based on the positioning error, the positioning error is held, the actual positioning observation result is input, and the positioning observation result is compared with the positioning error calculated by the positioning error simulation. A storage medium for storing a program for causing a computer to execute a positioning error evaluation method characterized by evaluating the following. Numerical simulation data representing the state of the atmosphere is input, stored in a memory, a gaze delay amount is calculated by an arithmetic unit based on the numerical simulation data, and the gaze delay amount is held. A positioning error is calculated by the positioning error simulation based on the positioning error, the positioning error is held, the actual positioning observation result is input, and the positioning error included in the positioning observation result is compared with the result calculated by the positioning error simulation. A program for causing a computer to execute a positioning error evaluation method characterized by evaluating.

Images