(57)【要約】【課題】測位誤差シミュレーションシステム、測位誤差評価方法、測位誤差シミュレーションを行うためのプログラム、及びそのプログラムを格納する記憶媒体に関する。大気の状態をほぼ正確に表している数値シミュレーションデータに基づいてＧＰＳ衛星方向の正確な視線遅延量を算出し、その視線遅延量により測位誤差シミュレーションをすることを目的とする。【解決手段】大気の状態を表す数値シミュレーションデータの入力手段と、数値シミュレーションデータに基づいて視線遅延量を計算する視線遅延量計算手段と、該視線遅延量計算手段で算出した視線遅延量に基づいて測位誤差を求める測位誤差計算手段とを備え、測位誤差をシミュレーションにより求める構成をもつ。【選択図】  図１
2004150901.tif 000002
【特許請求の範囲】
【請求項１】
大気の状態を表す数値シミュレーションデータの入力手段と、数値シミュレーションデータに基づいて視線遅延量を計算する視線遅延量計算手段と、該視線遅延量計算手段で算出した視線遅延量に基づいて測位誤差を求める測位誤差計算手段とを備え、測位誤差をシミュレーションにより求めることを特徴とする測位誤差シミュレーションシステム。
【請求項２】
数値シミュレーションデータは非静力学数値予報モデルに基づいて得られたものであることを特徴とする請求項１に記載の測位誤差シミュレーションシステム。
【請求項３】
実際の観測データに基づいて測定した測位観測結果の入力手段と、該測位観測結果と測位誤差シミュレーションで求めた測位誤差を比較し、実際の観測データにより得た測位観測結果の測位誤差を評価するための評価手段を備えることを特徴とする請求項１に記載の測位誤差シミュレーションシステム。
【請求項４】
データベース作成手段を備え、入力した数値シミュレーションデータの数値シミュレーションデータデータベースと、数値シミュレーションデータの作成条件を含む数値シミュレーションデータ情報を保持することを特徴とする請求項１に記載の測位誤差シミュレーションシステム。
【請求項５】
測位誤差シミュレーションで得た測位誤差を、その生成に使用した数値シミュレーションデータのモデルおよびマッピング関数で表したモデルを含む測位誤差生成情報とともにデータベースとして保持する測位誤差データベースを備えることを特徴とする請求項１に記載の測位誤差シミュレーションシステム。
【請求項６】
入力した実際の測位観測結果を測位観測結果データベースとし、測位観測結果情報とともに保持することを特徴とする請求項１に記載の測位誤差シミュレーションシステム。
【請求項７】
大気の状態を表す数値シミュレーションデータを入力し、メモリに保持し、該数値シミュレーションデータに基づいて演算手段により視線遅延量を算出し、該視線遅延量を保持し、演算手段により該視線遅延量に基づいて測位誤差シミュレーションにより測位誤差を算出し、該測位誤差を保持し、実際の測位観測結果を入力し、該測位観測結果を測位誤差シミュレーションにより算出した結果と比較することにより測位観測結果に含まれる測位誤差を評価することを特徴とする測位誤差評価方法。
【請求項８】
数値シミュレーションデータは非静力学数値予報モデルに基づいて生成したものであることを特徴とする請求項７に記載の測位誤差評価方法。
【請求項９】
入力した数値シミュレーションデータを数値シミュレーションデータの作成条件を含む数値シミュレーションデータ情報とともに数値シミュレーションデータデータベースとすることを特徴とする請求項７に記載の測位誤差評価方法。
【請求項１０】
測位誤差シミュレーションにおいて使用した数値シミュレーションデータのモデルおよびマッピング関数で表したモデルを含む測位誤差生成情報とともに該測位誤差をデータベースとして保持することを特徴とする請求項７に記載の測位誤差評価方法。
【請求項１１】
入力した実際の測位観測結果を測位観測結果情報とともに測位観測データデータベースとして保持することを特徴とする請求項７に記載の測位誤差評価方法。
【請求項１２】
大気の状態を表す数値シミュレーションデータを入力し、メモリに保持し、該数値シミュレーションデータに基づいて演算手段により視線遅延量を算出し、該視線遅延量を保持し、演算手段により該視線遅延量に基づいて測位誤差シミュレーションにより測位誤差を算出し、該測位誤差を保持し、実際の測位観測結果を入力し、該測位観測結果を測位誤差シミュレーションにより算出した該測位誤差と比較することにより測位観測結果を評価することを特徴とする測位誤差評価方法をコンピュータにより実行するプログラムを格納する記憶媒体。
【請求項１３】
大気の状態を表す数値シミュレーションデータを入力し、メモリに保持し、該数値シミュレーションデータに基づいて演算手段により視線遅延量を算出し、該視線遅延量を保持し、演算手段により該視線遅延量に基づいて測位誤差シミュレーションにより測位誤差を算出し、該測位誤差を保持し、実際の測位観測結果を入力し、該測位誤差シミュレーションにより算出した結果と比較することにより測位観測結果に含まれる測位誤差を評価することを特徴とする測位誤差評価方法をコンピュータにより実行するプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は測位誤差シミュレーションシステム、測位誤差評価方法、測位誤差シミュレーションを行うためのプログラム、及びそのプログラムを格納する記憶媒体に関するもので、特に、ＶＬＳＩもしくはＧＰＳによる測位誤差シミュレーションをするものである。大気の状態をほぼ正確に表している数値シミュレーションデータに基づいてＧＰＳ衛星方向の正確な視線遅延量を算出し、その視線遅延量により測位誤差シミュレーションをすることができ、正確な測位誤差を求めることができる。数値シミュレーションデータは数値シミュレーションにより大気の状態を数値データとして求めたものである。また、その測位誤差シミュレーションにより求めた測位誤差をもとに実際の観測データに基づいて得た測位観測結果と比較し、実際の測位観測結果の測位誤差を容易に評価することができるようにした。
【０００２】
【従来の技術】
従来、位置測定は人工衛星もしくはクエーサー等から発信される電波を使用し、ＶＬＳＩであれば二カ所の基地における電波の観測時間の差にもとづいて位置を測定している。ＧＰＳであれば人工衛星から発信される電波の到達時間から衛星までの距離を複数の人工衛星について推定し位置を測定する。
【０００３】
人工衛星もしくはクエーサーから発信される電波は、大気中で屈折率により伝搬速度が減少することと、伝搬経路が曲率をもち、真空中での伝搬経路より見かけ上長くなり大気伝搬遅延を生じる。図１３（ａ）はＧＰＳにより測位する場合に、ＧＰＳから放射される電波が大気中を伝搬する際に生じる伝搬時間の遅延（以下、視線遅延量と称する）の説明図である。ＧはＧＰＳ衛星、クエーサーなどの電波源から地上のアンテナまでの直線距離である。Ｓは大気屈折により生じる実際の経路を表す。１３（ｂ）は天頂遅延量とマッピング関数に基づいて視線遅延量を推定する方法の説明図である（後述する）。
【０００４】
図１４は、従来の測位観測システムを示す。図１４において、４０１は観測データであって、ＧＰＳ衛星から測定点Ｐまでの電波伝搬時間、ＧＰＳ衛星方位、仰角等である。４０２は測位計算システムであって、観測データに基づいて測位観測により測位誤差を計算するものである。４１０は測位計算手段である。４１２は天頂遅延量推定の処理であり、観測データ（電波伝搬時間）とマッピング関数により天頂遅延量を推定する処理である。４１３は測位誤差算出の処理であって、視線遅延量（天頂遅延量）に基づいて測位誤差を計算するものである。４２０はマッピング関数である。４２１は測位観測結果の出力を表す。
【０００５】
従来から行われている実施の観測データにより行われる測位誤差の観測は、観測された伝搬時間を基に天頂遅延量を求め、その視線遅延量に基づいて測位誤差を計算することにより行なわれていた。
【０００６】
【発明が解決しようとする課題】
上記のように、従来の測位誤差の観測は、大気中を伝搬する実際の観測データに基づいて行われ、大気伝搬における屈折率による伝搬速度の遅れ、伝搬経路に生じる曲率による誤差などのために、大気の影響により測定する毎に異なって観測される。その変動が大気の状態によるものであるのか、あるいは、地核変動等のために実際に観測位置が移動したものであるのかを判定することは困難であった。
【０００７】
本発明は、大気の状態を正確に表している数値シミュレーションデータに基づいて正確な視線遅延量を計算し、正確な測位誤差を求めることのできる測位誤差シミュレーションシステムを提供し、実際の測位観測結果と比較することにより実際の観測データで観測された測位誤差を容易に評価することができるようにすることを目的とする。
【０００８】
【課題を解決するための手段】
気象学における数値シミュレーションデータは、大気空間を３次元の立方体の格子空間に分割し、気象観測データをもとに各格子の大気の状態を計算したものである。格子間隔を小さくし（数百メートルから数キロメートル）、非静力学モデルで解析したものは、大気の状態をかなり正確に表現している。本発明は、そのような気象学の数値シミュレーションデータに基づいてＧＰＳ衛星から観測点にいたる電波の視線遅延量を正確に計算し、その視線遅延量に基づいて測位誤差シミュレーションを行うことにより、正確に測位誤差を計算できるようにした。
【０００９】
本発明の実施の形態を説明するのに先立ち、大気屈折の基本原理について説明する。大気屈折による大気遅延は次の２つの物理的効果に起因する。１つは、マイクロ波が誘電媒質である中性大気中を通過するために真空中よりも減速される効果である。２つ目は、伝搬経路が曲率（  ｒａｙｂｅｎｄｉｎｇ）を持つため、直線より経路が実際に伸びる効果である。このとき、大気の密度は下層ほど大きいので、上方に凸の伝搬経路となる。これらの効果による遅延ΔＬをＥＰＤ（Ｅｘｃｅｓｓ  Ｐａｔｈ  Ｄｅｌａｙ）と呼び、次式で表される。
【００１０】
【数１】
2004150901.tif 000003
【００１１】
ここで、ＧはＧＰＳ衛星・クエーサーなどの電波源から地上のアンテナまでの直線距離、ｎ（ｓ）は大気屈折によって曲率が生じた伝搬経路Ｌ上の点ｓにおける屈折率である。このとき、Ｓの長さはＬ上の微少部分ｄｓを積分して得られる。右辺第１項が減速の効果であり、実際の電波はＬ上を中性大気の屈折率分だけ遅い伝搬速度で伝搬する。ところが、ＶＬＳＩやＧＰＳの解析の上では同じ経路上を真空中の速度で伝搬したとみなすため、双方の速度差に伝搬時間を乗じた分だけ見かけ上伝搬経路が伸びたことになる。残りの右辺第２項が伝搬経路の曲率の効果となり、当然の事ながらこの効果による遅延は天頂方向で０となる。
【００１２】
ここで、湿潤大気の屈折率の計算は以下のような式で得られる（  ここでは途中の過程を一部省略した）  。
【００１３】
【数２】
2004150901.tif 000004
【００１４】
ここで、式（２）のＰは全大気圧（ｈＰａ）、Ｐｖ及びＴはそれぞれ水蒸気分圧（ｈＰｗ）、及び絶対温度（Ｋ）である．各項の係数Ｋ１  、Ｋ２  室内実験より決定される定数であり、研究者によって若干異なる。一般に、式（２）の右辺第１項を“静水圧項”と呼ぶ。この項は、静水圧平衡にある大気成分の密度に比例する寄与であり、これを地上から大気上層まで積分すると静水圧遅延が得られる。式（２）の右辺第２項は“湿潤項”と呼ばれ、これを積分して得られる遅延を湿潤遅延と呼ぶ。
【００１５】
今回の非静力学モデルによるレイトレーシング計算では、まず非静力学モデルの各格子点に与えられた気温、気圧、水蒸気分圧の各値から式（２）を用いて格子点毎の屈折率を計算する。次に、地上のある観測点から任意の仰角と方位で射出した電波が衛星まで到達する経路を、計算した屈折率からＳｎｅｌｌの法則に従って逐次屈折角を計算する。この計算そのものがレイトレーシングである。レイトレーシングで計算された視線方向の遅延量と、従来の大気モデル（これがマッピング関数と呼ばれるもの）  とを比較することで、従来のモデルの評価ができる。次に本発明の測位誤差シミュレーションについて説明する。まず従来の大気推定法であるマッピング関数について簡単に説明し、その後でレイトレーシングした結果による測位誤差シミュレーションを説明する。
【００１６】
従来の大気遅延推定、即ち、マッピング関数ＧＰＳやＶＬＳＩの解析では、２つの観測局間での電波の位相差データに含まれる観測局間の相対位置、電離層遅延、大気遅延、観測局ないしは人工衛星内の時計の誤差等の寄与を未知パラメータとして膨大な数の観測方程式を立て、各パラメータを最小二乗法で推定する手法を取る。大気屈折に関していえば、厳密には気温・気圧・水蒸気分圧の時間的・空間的変動が存在するため、視線方向の大気遅延は時間・空間の双方の関数となる。しかし、これらをすべて未知数とするとパラメータが多すぎて観測方程式を解くことが不可能になる。そこで、ある一定の時間内において大気が水平方向に一様な構造を維持すると仮定し、その間は大気遅延が仰角のみに依存するとするモデル化を行えば未知数を大幅に減らすことができる。またモデル化により、デジタル化された観測データを数値計算する上で大気遅延の扱いが容易となる。この思想で導入されたのがマッピング関数と呼ばれるモデルである。マッピング関数を使うと、ある仰角θの全大気遅延ΔＬは下記のように簡単に表すことができる。
【００１７】
【数３】
2004150901.tif 000005
【００１８】
ここで、ΔＬｚ  ｈ  とΔＬｚ  ｗ  は、それぞれ天頂方向の静水圧遅延と湿潤遅延である。これらに乗じられているＭｈ  （θ）とＭｗ  （θ）が、それぞれ静水圧遅延用と湿潤遅延用のマッピング関数である。
【００１９】
後でも触れるが、マッピング関数は、観測方程式を線形化した際に作られる計画行列の中の大気遅延に関する偏微分係数そのものとなる。ただし、ＶＬＳＩでは２地点間でのマッピング関数の差（後述する「一重位相差観測値」に含まれる“大気遅延量差”）、ＧＰＳではさらに異なる衛星について求めた２地点間でのマッピング関数の差どうしで再び差を取った値（後述する「二重位相差観測値」に含まれる“大気遅延量差の差”）となる。
【００２０】
マッピング関数は、基本的には仰角θの正弦の逆数、すなわち（１／ｓｉｎ（θ））となる。ただし、１０度以下の低仰角の観測を行う場合には、マイクロ波が大気中でより長い経路をたどるため、伝搬経路の曲率がきつくなる。その結果、曲率の効果による遅延が指数関数的に増大する。しかも、地球の曲率の影響も無視できなくなる。そこで、単純な「仰角の正弦の逆数」ではなく、下記の式（４）に示すように分母を連分数で拡張したマッピング関数が基本形として考えられている。
【００２１】
【数４】
2004150901.tif 000006
【００２２】
ここで、各項の係数（ａ、ｂ、ｃ）は定数である。これらの係数は、球対称大気モデルから得られる屈折率プロファイルからレイトレーシングを用いて求めた大気遅延とマッピング関数とを比較し、双方が一致するように最小二乗法を用いて決められる。実際のマッピング関数の開発では、係数を常数としたり、地上気象データの気圧・温度などで変化する変数にしたり、また最後のｂ／（ｓｉｎ  θ  ＋ｃ）の部分をｔａｎθなど他の関数型に置き換えたりなど、研究者によってより高精度のマッピング関数を求めて様々な改良がなされている。こうして開発されたマッピング関数と、波線追跡法から、それぞれ計算された大気遅延量の差の大小が、マッピング関数の適用範囲、すなわち推定精度を示すことになる。
【００２３】
具体的な手順としては、ＧＰＳやＶＬＳＩ解析で使用される観測方程式に、レイトレーシングした遅延量を入力して観測方程式を解くことにより、大気遅延によって生じる測位誤差シミュレーションが可能になる。
【００２４】
具体的には次のように計算する。まず、ＧＰＳやＶＬＳＩの計算では２地点間の位相差を取ることから、遅延量についても観測点間での差を計算する。観測点Ｒ１  とＲ２  に共通の衛星から電波が平行に到達すると仮定すると、観測エポックでｉ番目の２地点間の遅延量差は、
【００２５】
【数５】
2004150901.tif 000007
【００２６】
ここでΔＬＳ  （Ｒ１  ）ｉ  とΔＬＳ  （Ｒ２  ）ｉ  は、それぞれ観測点Ｒ１  とＲ２  における任意の方位・仰角での視線遅延量である。添え字のＳとｉはそれぞれ衛星とエポックの番号を示す（  エポックとは観測の区切りである）  。ここで、観測点Ｒ１  が固定点とすると、遅延量差ΔＬＳ  （Ｒ２  −Ｒ１  ）ｉ  は２地点間の遅延量差であり、通常のＧＰＳ解析で言うところの「一重位相差観測値」に含まれる大気遅延量差に相当する。
【００２７】
ここで、多数のエポック毎の遅延量差を観測データとして、この遅延量差に最も適合するベクトルを最小二乗法で計算することができる。すなわち、この関係を線形化すると、
【００２８】
【数６】
2004150901.tif 000008
【００２９】
である。式（６）で、行列ｅは誤差ベクトル、行列Ｙが遅延量残差から得られる、エポック１番目からｍ番目までの観測量のベクトルとおけて、
【００３０】
【数７】
2004150901.tif 000009
【００３１】
と表される。また、行列Ａは、各要素の偏微分係数からなる計画行列であり、簡単のため推定するパラメータを観測点位置の誤差ベクトル（後述する見かけの”変位量ベクトル”）  と天頂方向の大気遅延量とすると、
【００３２】
【数８】
2004150901.tif 000010
【００３３】
となる。ここで、
αｉ  、βｉ  およびγｉ  （  ｉ＝１、２、・・・・ｍ）  などは、遅延量差のｘ、ｙ、ｚの各成分についての方向余弦である。また、Δｍｆ（θ）  （Ｒ２  −Ｒ１  ）は遅延量についての偏微分係数である（  ここではθを仰角とするマッピング関数の２地点間での差）  。さらに、Ｄが、この方程式であてはめを行った結果得られるベクトルであり、
【００３４】
【数９】
2004150901.tif 000011
【００３５】
と表現できる。ここで、ｘ、ｙ、ｚはベクトルの各成分であり、ｘ軸の正の方向に北、ｙ軸の正の方向に東、ｚ軸の正の方向に鉛直上向きとする。また、ＺＤは天頂方向の遅延量である。
【００３６】
このベクトルＤは、観測点Ｒ２  の位置を真の位置からずらすようなベクトルであり、これを見かけの変位量ベクトル、すなわち大気遅延によって生じる測位誤差と見なすことができる。この見かけの変位量ベクトルの計算が測位誤差シミュレーションであり、様々な大気状態の元でこのベクトルの大きさや方向を調べることで測位誤差の振る舞いを調べることが可能になる。
【００３７】
実際に、最小二乗法を用いて遅延量差から見かけの変位量ベクトルＤを求めることを考えると、この場合の観測方程式は、
【００３８】
【数１０】
2004150901.tif 000012
【００３９】
のように書ける。
【００４０】
ここで得られる行列Ｄは、例えばＶＬＳＩ解析における大気遅延によって生じる見かけの変位量に相当する。さらに実際のＧＰＳ解析では、２地点の観測局に置かれたＧＰＳ受信機内の時計誤差を相殺するために、異なる衛星についての位相差どうしをさらに引き算する。これを「二重位相差観測値」という。そこで、測位誤差シミュレーションを行う場合についても同様に、遅延量差の差を計算して見かけの変位量ベクトルＤを計算する。２つの異なる衛星Ｓ１  とＳ２  について遅延量の差の差（二重遅延量差）  から行列Ｄを求めるときの行列Ｙは次のように書ける。
【００４１】
【数１１】
2004150901.tif 000013
【００４２】
【数１２】
2004150901.tif 000014
【００４３】
であり、計画行列Ａは、
【００４４】
【数１３】
2004150901.tif 000015
【００４５】
と表される。より現実に近い測位誤差シミュレーションを行う場合は、４つないしそれ以上の衛星を想定して最小二乗法を実施する。すなわち、ｎ個の衛星について計算する場合は、行列ＹとＡの列の数が、一重差の場合はｎ×ｍ列となり、二重差の場合は（ｎ−１）×ｍ列となる。
【００４６】
また、先の式（８）で、仰角θの他に方位Ｒに依存した遅延勾配を考慮したマッピング関数ｍｆ（θ，  α）  （遅延勾配モデルとも言う）を大気遅延推定に用いたときのシミュレーションももちろん可能である。このときは、偏微分係数となるｍｆ（θ，  α）  は、大気の非一様性を単純な平面勾配で近似して構築され、ＭａｃＭｉｌｌａｎ（１９９６）を例に挙げると以下の式で表される。
【００４７】
【数１４】
2004150901.tif 000016
【００４８】
右辺第１項が従来の球対称大気モデルのマッピング関数、右辺第２項が大気の非一様性の効果を示す。この第２項を求める手順としては、まず平面近似された大気勾配の大きさと向きをベクトルＧ（勾配ベクトル）  で表す。上式で、Ｇの南北成分がＧＮ  、東西成分がＧＥ  である。次に電波の到来方位αの単位ベクトルを考慮して大気遅延の水平勾配の強度と方位を表す量を求める。言い換えれば、この操作ではＧと電波到来方向のベクトルとの内積をとることになる。この第２項で乗じられているｃｏｔ（θ）  は、仰角が減少するにしたがって水平方向の伝搬経路が長くなる結果、水平勾配の強度が強調される効果を考慮した項、及びｍ０  （θ）  は以上のすべての効果を天頂遅延に投影するための項である。ｍ０  （θ）  は従来の球対称マッピング関数を用いる。ここで、本来は第２項のｃｏｔ（θ）  に含まれるθの大きさは、大気の非一様性を考慮したときの仰角θ０  を用いるべきである。しかしながら、双方の差（θ’−θ）  は仰角５度でも０．０３度程度に過ぎないため、近似的にθが使用されている。
【００４９】
ここでまとめると、遅延勾配モデルを評価対象として測位誤差シミュレーションを行う場合、天頂遅延量ＺＤの他に勾配ベクトルＧの南北成分ＧＮ  ，東西成分ＧＥ  、の２つが推定パラメータとして加わることになる。
【００５０】
こうしてシミュレートされた測位誤差は、次のようにして実際の測位解の評価に使用できる。つまり、実際のＧＰＳやＶＬＳＩ測位がなされた同時期、同地域の非静力学モデルから計算した測位誤差と、実際の測位結果のばらつきと比べることで、測位結果の評価が可能になる。例えば、突発的な地殻変動に見えるような変位がＧＰＳによって解析された場合に、測位誤差シミュレーションで同様の変位が計算されれば、実際の地殻変動ではないと見なされるし、逆であれば地殻変動の可能性が高いと判断できる。また、実際の測位解の振る舞いと、シミュレートされた測位誤差との比較を様々な気象条件下で行うことによって、現行の大気モデルの改良に必要な基礎データを取得することができる。
【００５１】
図１は本発明の概念の説明図である。図１において、Ａは気象学の数値シミュレーションデータによる大気遅延量（視線遅延量）、測位誤差の評価システムを表す。Ｂは既存の解析システムであって、ＧＰＳもしくはＶＬＳＩによる測位観測システムである。本発明のシステムＡは、気象学に用いる数値シミュレーションデータ５０１をもとに数値計算５０２により、視線遅延量および測位誤差シミュレーションを行う。５０３は数値シミュレーションデータをもとにレイトレーシングによる視線遅延量計算の処理を表す。５０４は数値シミュレーションデータに基づいて５０３の処理で求めた視線遅延量により正確な測位誤差をシミュレーションにより求め、その測位誤差と実際の観測データにより測位観測結果と比較することにより実際の測位結果を評価する測位誤差シミュレーションの処理を表す。
【００５２】
既存の解析システムＢは、図１４で説明したものと同様である。球対称モデル、あるいは遅延勾配モデル等のモデルに基づく大気遅延量推定の処理４２３は、それぞれのモデルのマッピング関数により大気遅延量を推定する処理を表す（球対称モデルおよび遅延勾配モデルによる大気遅延量推定は図１５参照）。
【００５３】
図２は、本発明のシミュレーションによる測位誤差算出方法の概要を示す図である。図２において、６１０は、本発明で導入可能なベクトルＹの種別を表す。６０１は観測方程式を表す。６０２は観測量のベクトルＹを求める処理を表し、視線方向の遅延量の計算結果ΔＬ等である。６０３は未知パラメータのベクトルＤを表し、視線遅延量によって生じる見かけの変位（測位誤差）ベクトル、天頂遅延量ＺＤ等を要素とするものである。
【００５４】
６０４は計画行列Ａを表し、各パラメータの偏微分係数が要素であるものである。見かけの変位（シミュレーションにおける測位誤差）ベクトルは各成分の方向余弦である。天頂遅延量はＺＤであり、勾配遅延パラメータＧＮ  、ＧＥ  についてはマッピング関数で表す。６１２は最小二乗法によりベクトルＤを推定する処理を表す。
【００５５】
上記の観測方程式（１０）は既存の解析システムＢにおいても同じである。しかし、従来の測位誤差シミュレーションと本発明の測位誤差シミュレーションではそこで扱うパラメータが以下のように、行列Ｙ（ベクトルＹ）、行列Ｄ（ベクトルＤ）において相違する。
【００５６】
実際の測位において既存の解析システムＢにおける行列Ｙ、行列Ｄは次のようなものである。行列Ｙにおいて、エポック１番目からｍ番目まで、それぞれのエポックでのＧＰＳ衛星の位置（軌道情報ＸＳ  、ＹＳ  、ＺＳ  ）から地上の観測点（ＸＲ  、Ｙ  Ｒ、ＺＲ  ）  に到達した電波の伝搬時間に光速を乗じて計算した距離である。この量には、地上の観測点と衛星との間の真の距離情報、大気遅延量、衛星の時計誤差、観測点に設置した受信機の時計誤差などが含まれるため、真の距離と区別して疑似距離と呼ぶ。行列Ａは各推定要素の偏微分係数からなる計画行列である。
【００５７】
行列Ｄの未知パラメータにおいて、観測点位置（ｘ、ｙ、ｚ）  、及び大気による天頂方向の天頂遅延量（ＺＤ）は次のようなものである。既存の解析システムＢでは、その推定されるＺＤはあくまで、最小二乗法によってつじつまが合うように数学的に求められる値であり、実際の大気の状態がどのようになっているかとは無関係なものである。また、既存の解析システムＢでは大気中の水蒸気の変動が激しいときなどは、このＺＤの推定精度は落ちるものである。
【００５８】
本発明における測位誤差シミュレーションにおける行列Ｙ、行列Ｄは次のようなものである。行列Ｙは、エポック１番目からｍ番目まで、それぞれのエポックでのＧＰＳ衛星の位置（方位・仰角）をソフトウェアに入力し、数値シミュレーションデータを用いてレイトレーシングで計算した視線方向の大気遅延量である。したがって、本発明の測位誤差シミュレーションでは数値シミュレーションデータから求めた正確な大気状態に基づく伝搬遅延量が計算される。
【００５９】
本発明の測位誤差シミュレーションにおける行列Ｄの測位誤差ベクトル、及び大気による天頂方向の天頂遅延量（ＺＤ）は以下のようなものである。測位誤差ベクトルは、未知パラメータである観測点の真の位置（ＸＲ  、ＹＲ  、ＺＲ  ）から大気遅延によってずれた位置（ＸＡ  、ＹＡ  、ＺＡ  ）までのベクトルである。すなわち、本発明の測位誤差シミュレーションにおいては、数値シミュレーションデータから求めた大気のみによって生じる伝搬遅延量が計算され、入力されているので、ここでの測位誤差ベクトルや天頂遅延量ＺＤは、明らかに大気によって生じた量と見なすことができる。したがって、本発明の測位誤差シミュレーションでは、推定される測位誤差や天頂遅延量は真の値に近いものである。また、本発明で計算した測位を使用すれば、入力する数値シミュレーションデータそのものの精度も向上し、より現実の大気の状態を推定することが可能になる。
【００６０】
さらに、実際の測位観測結果と本発明の測位誤差シミュレーションで推定された測位誤差とを比較して、双方の変動のパターンがほぼ同じであれば、実際の測位観測結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測結果の変動が地殻変動によるものと判定できる。
【００６１】
上記のように、本発明によれば、気象学の数値シミュレーションデータに基づいて正確に視線遅延量を計算し、その視線遅延量に基づいて測位誤差シミュレーションを行うことができる。そのため正確な測位を容易に把握することができるようになる。また、実際の観測で得た測位観測結果との相関や、そのデータを時系列データとして配列してその傾向を見ることが容易になり、測位観測結果のバラツキ等の原因を容易に判断できるようになる。また、実際の観測データをもとに測位計算する解析システムで採用したマッピング関数の精度、測位誤差の評価等も容易に行うことができる。
【００６２】
【発明の実施の形態】
図３は本発明の測位誤差シミュレーションシステムの構成図である。図３において、１はＣＰＵであり、測位誤差シミュレーションの演算を行うものである。２はメモリであり、測位誤差演算において必要なデータを記憶するものである。３は入出力インターフェイスである。４はデータ入力装置である、本発明で使用する気象学の数値シミュレーションデータ等を入力するものであり、６はデータベース作成手段であって、本発明で使用する数値シミュレーションデータ、本発明の測位誤差シミュレーションで求めた測位誤差等をデータベース化するものである。７は測位誤差シミュレーション手段であり、視線遅延量計算手段１１、測位誤差計算手段１２、評価手段１３を備えるものである。８はデータ保存装置であり、数値シミュレーションデータデータベース２１、マッピング関数保持手段２２、測位誤差シミュレーション結果データベース２３、測位観測結果データベース２４、各種データ保持手段２５を備えるものである。
【００６３】
図４は本発明の測位誤差シミュレーションシステムの装置構成の詳細である。図４の構成において、視線遅延量計算手段１１は前述の式（１）、（２）、（３）と（４）の数式の演算手段を備え、レイトレーシング法によりＧＰＳ方向の視線遅延量および天頂遅延量を求めるものである。測位誤差計算手段１２は前述の行列Ｙ、行列Ａ、行列Ｄにより最小二乗法により観測方程式の演算を行うものである。測位観測システム５１は既存のもので、図１の解析システムＢと同じである。
【００６４】
図４の構成において、数値シミュレーションデータ選択手段３４は、他システム３１から数値シミュレーションデータを入力する。数値シミュレーションデータは、例えば、非静力学モデルで計算されたものであり、大気の状態を正確に表しているもので、例えば、インターネット等のネットワークを介して入力する。あるいは、記憶媒体に格納されている状態で入力することもできる。もしくは、本発明のシステムに備える数値シミュレーションデータデータベース２１から入力する。入力された数値シミュレーションデータはメモリ領域３７に保持される。視線遅延量計算手段１１は数値シミュレーションデータとレイトレーシング法によりＧＰＳ方向の視線遅延量および天頂遅延量を計算する。計算された視線遅延量は視線遅延量保持手段４１に保持される。測位誤差計算手段１２は、数値シミュレーションデータをもとに求めた視線遅延量をもとに、行列Ｙ、行列Ａ、行列Ｄにより観測方程式の演算を行い、測位誤差を求める。そのようにして求めた測位誤差は測位誤差シミュレーション結果保持手段４２に保持される。評価手段１３は測位誤差シミュレーション結果の測位誤差と、既存の測位観測システム５１で得た測位誤差とを入力し、比較する。その統計的処理を行い、既存の測位観測誤差のバラツキ等の原因を視覚的に評価できるようにする。出力手段５４は評価手段１３で求めた評価結果をディスプレイに表示する、もしくはプリンタに印刷出力する。あるいは、記憶媒体に保存する。
【００６５】
入力した数値シミュレーションデータはデータベース作成手段３８により、数値シミュレーションデータ情報とともに数値シミュレーションデータデータベース２１および数値シミュレーションデータ情報保持手段３５に保存する。本発明の測位誤差シミュレーションでもとめた測位誤差は測位誤差シミュレーション結果データベース２３に保存する。また、既存の測位観測システム５１で得た測位誤差（測位観測データ）は測位観測結果データベース２４に保存することもできる。
【００６６】
図５は本発明の視線遅延量算出のための動作説明図である。本発明の視線遅延量算出は、数値シミュレーションデータとＧＰＳ衛星の位置に基づいて前述の式（１）、（２）、（３）の計算を行うものである。レイトレーシング法により視線方向（ＧＰＳ方向）の遅延量を求める。また天頂遅延量をもとめる。天頂遅延量は、数値シミュレーションデータをもとに計算する方法と視線遅延量をもとにマッピング関数を使用して求める等さまざまな方法で計算できる。
【００６７】
図５の構成の動作を説明する。視線遅延量計算手段１１は、数値シミュレーションデータ入力処理１１１により数値シミュレーションデータを入力する。数値シミュレーションデータとしては、気温、気圧、水蒸気量（露点温度）であり、各データは位置を示す格子点、年月日、時刻をもつ。入力された数値シミュレーションデータはメモリ領域１２０に保持される。屈折率算出手段１１２において、静水圧項計算手段１１３により静水圧項を計算する。湿潤項計算手段１１４は数値シミュレーションデータの気温、気圧、水蒸気量（露点温度）に基づいて、式（２）により湿潤項を計算する。静水圧項遅延量算出手段１１５と湿潤遅延量算出手段１１６により式（１）の第一項の計算をし、静水圧遅延量と湿潤遅延量をメモリ領域１２０に保持する。また、（Ｓ−Ｇ）計算手段１１７は数値シミュレーションデータおよびＧＰＳ衛星位置から式（１）の（Ｓ−Ｇ）の計算をし、メモリ領域１２０に保持する。静水圧遅延量と湿潤遅延量および（Ｓ−Ｇ）により視線方向遅延量ΔＬを求める。またレイトレーシング法により天頂遅延量を求める。天頂遅延量は視線遅延量をもとにマッピング関数で求める方法もある。使用するマッピング関数はマッピング関数指定入力１３１により指定され、選択される。算出した視線遅延量は視線遅延量保持手段４１に保持する（天頂遅延量の保持データは図示を省略している）。
【００６８】
図６は本発明の測位誤差算出のための動作説明図である。本発明の測位誤差シミュレーションは数値シミュレーションデータにより計算した視線遅延量をもとに、ＶＬＳＩでは前述の式（５）、（６）、（７）、（８）、（９）、（１０）により測位誤差シミュレーションを行う。ＧＰＳでは、前述の式（５）、（６）、（９）、（１０）、（１１）、（１２）、（１３）により測位誤差シミュレーションを行う。
【００６９】
一重位相差により測位誤差シミュレーションを行うか、二重位相差による測位誤差シミュレーションを行うか決め、また使用する大気モデルを決めて、測位誤差計算手段１２に数値シミュレーションデータで求めた各エポックの衛星方向の視線遅延量を入力し、行列Ｙを計算する（ベクトルＹ演算１５１の処理）。求めた行列Ｙはメモリ領域１６０に保持する。また、計画行列Ａを生成し（行列Ａ計算１５３の処理）、メモリ領域１６０に保持する。行列Ｄ（位置座標ｘ、ｙ、ｚ、天頂遅延量ＺＤが未知パラメータ）を設定し（未知パラメータＤ設定の処理１５２）、メモリ領域１６０に保持する。最小二乗法により誤差ベクトルｅが最小になるように、観測方程式（１０）に従って行列Ｄ（  ｘ、ｙ、ｚ、ＺＤ）  を求める。このようにして求めたｘ、ｙ、ｚは測位誤差であり、ＺＤは推定天頂遅延量である。求めた測位誤差、推定天頂遅延量は測位誤差シミュレーション結果保持手段１８０に保持する。求めた測位誤差シミュレーションで求めた測位誤差はデータベース作成手段３８により測位誤差シミュレーション結果データベース２３に保存される。
【００７０】
測位誤差を評価する場合には、本発明に測位誤差シミュレーションで求めた測位誤差と実際の測位観測結果を評価手段１３に入力する。評価手段１３は測位誤差シミュレーションにより推定された測位誤差とを比較し、相違を視覚的にわかるように出力する。統計的に処理した双方の変動のパターンがほぼ同じであれば、実際の測位解析結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測システムで得た測位観測結果の変動が地殻変動によるものであると判定できる。
【００７１】
図６において、実際の測位観測システムから測位観測結果を入力する場合には、測位観測データ情報保持手段１７１を参照し、選択された測位観測データ情報選択手段１７２により測位観測データを選択して、測位観測結果データベースから入力する。
【００７２】
図７は、既存の解析システムで観測した測位観測結果をデータベースにするための構成である。測位観測結果入力手段２０１により測位観測システム２００から測位観測結果を入力する。データベース作成手段３８により、測位観測結果データベース情報を作成し、測位観測結果情報保持手段５２に保持するとともに、測位観測結果は測位観測結果データベース２４に保存する。
【００７３】
図８は本発明の視線遅延量算出のフローチャートである。図４、図５を参照し、図８の説明をする。数値シミュレーションデータ選択手段３４で選択された数値シミュレーションデータが入力され、数値シミュレーションデータを作成したモデル、数値シミュレーションデータの対象とする年月日、時刻、位置等を選択する（Ｓ１）。選択した数値シミュレーションデータを入力し、保持する（Ｓ２、Ｓ３）。選択した数値シミュレーションデータをデータベースに登録する場合には、データベース登録処理をし、データベースに登録する（Ｓ４、Ｓ５、Ｓ６、Ｓ７）。Ｓ５でデータベース登録の処理を行わないのであれば、数値シミュレーションデータから気圧、気温、水蒸気分布を求める（Ｓ８）。静水圧項の係数Ｋ１  、湿潤項の係数Ｋ２  を定める（Ｓ９）。静水圧項Ｋ１  （Ｐ／Ｔ）と湿潤項Ｋ２  （Ｐｕ／Ｔ２  ）を計算する（Ｓ１０）。屈折率に関する式（ｎ−１）×１０６  を計算する（Ｓ１１）。（Ｓ−Ｇ）の計算をする（Ｓ１１’）。静水圧遅延ΔＬｈ  と湿潤遅延ΔＬｗ  、それぞれの天頂遅延量ΔＬｈ  ｚ  、ΔＬｗ  ｚ  を計算する（Ｓ１２）。静水圧遅延ΔＬｈ  と湿潤遅延ΔＬｗ  、天頂遅延量ΔＬｈ  ｚ  、ΔＬｗ  ｚ  を保持する（Ｓ１３）。
【００７４】
図９は本発明の測位誤差算出のフローチャートである。図９は一重位相差の場合を示す。図９の説明で図６を参照する。まず、ΔＬｈ  とΔＬｗ  を入力する（Ｓ１）。各イベント（各位置、視線方向）で、式（５）、もしくは式（１２）を計算する。ベクトルＹ（行列Ｙ）を計算する（Ｓ３）。計画行列Ａを求める（Ｓ４）。ベクトルＤを設定する（Ｓ５）。最小二乗法により誤差ベクトルｅを最小にするベクトルＤのパラメータ値を計算する（Ｓ６）。計算したベクトルＤの行列要素を測位誤差として保存する（Ｓ７）。
【００７５】
図１０は本発明の評価手段の視線遅延量の評価のフローチャートである。評価手段は、本発明の測位誤差シミュレーションシステムで求めた視線遅延量と、実際の測位観測システムでマッピング関数をもとに推定した視線遅延量を入力する（Ｓ１）。実際の測位観測システムで推定し視線遅延量と測位誤差シミュレーションで算出した視線遅延量を比較する（Ｓ２）。計算された誤差の可視化、統計処理を行う（Ｓ３）。可視化されたデータや統計処理結果を保存する（Ｓ４）。
【００７６】
図１１は本発明の評価手段の測位誤差の評価のフローチャートである。測位誤差シミュレーションで求めた測位誤差を入力する（Ｓ１）。測位誤差シミュレーションで求めた結果と実際の測位観測システムの観測結果を比較し、双方の相関や時系列のデータを可視化する（Ｓ２）。可視化されたデータや時系列データを保存する（Ｓ３）。
【００７７】
図１２は、本発明の各種データの構成を示す。（ａ）は数値シミュレーションデータの例を示す。（ｂ）は本発明における測位誤差シミュレーションにおける視線遅延量のデータ構成の例を示す。各ＧＰＳ衛星についてもつものである。（ｃ）は本発明における測位誤差シミュレーションによる測位誤差の例を示す。各ＧＰＳ衛星についてもつものである。（ｄ）は実際の測位観測システムにおける測位観測データの例を示す。
【００７８】
【発明の効果】
本発明の測位誤差シミュレーションにおいては、数値シミュレーションデータから求めた大気のみによって生じる伝搬遅延量を計算することができる。したがって、本発明の測位誤差シミュレーションでは、大気のみにより推定される測位誤差や大気に依存しない天頂遅延量を得ることができる。また、本発明では計算した視線遅延量を用いることにより、高精度のマッピング関数を考案できる。
【００７９】
さらに、実際の測位観測システムにより得られた測位観測結果と本発明の測位誤差シミュレーションで推定された測位誤差とを比較したとき、双方の変動のパターンがほぼ同じであれば、実際の測位観測結果に大気による誤差が顕著に含まれることを意味し、パターンが似ていなければ、実際の測位観測結果の変動が地殻変動によるものと判定できる。
【００８０】
このように、本発明によれば、正確な測位を容易に把握することができるようになる。また、実際の観測で得た測位観測結果との相関や、そのデータを時系列データとして配列してその傾向を見ることにより、測位観測結果のバラツキ等の原因を容易に判断できるようになる。また、実際の観測データをもとに測位計算する解析システムで採用したマッピング関数の精度、測位誤差評価等も容易に行うことができる。
【図面の簡単な説明】
【図１】本発明の概念の説明図である。
【図２】本発明の測位誤差算出方法の概要を示す図である。
【図３】本発明の測位誤差シミュレーションシステムの構成図である。
【図４】本発明の測位誤差シミュレーションシステムの装置構成の詳細を示す図である。
【図５】本発明の視線遅延量算出のための動作説明図である。
【図６】本発明の測位誤差算出のための動作説明図を示す。
【図７】既存の測位観測システムの測位観測結果についての本発明のデータベース作成の動作説明図を示す。
【図８】本発明の視線遅延量算出のフローチャートを示す図である。
【図９】本発明の測位誤差算出のフローチャートを示す図である。
【図１０】本発明の視線遅延量の評価のためのフローチャートを示す図である。
【図１１】本発明の測位誤差評価手段のフローチャートを示す図である。
【図１２】本発明の各種データの構成を示す図である。
【図１３】電波の大気伝搬遅延とマッピング関数の説明図を示す。
【図１４】既存の測位観測システムを示す図である。
【図１５】大気モデルの模式図を示す。
【符号の説明】
Ａ：大気遅延量・測位誤差評価システム
Ｂ：既存の解析システム
４０１：観測データ
４０２：測位計算システム
４１０：測位計算
４２３：大気遅延量推定
４２６：測位結果
５０１；気象学に用いる数値シミュレーションデータ
５０２：数値計算
５０３：レイトレーシングによる視線遅延量計算
５０４：実際のＧＰＳやＶＬＳＩ解析を考慮した測位誤差シミュレーション
【図１】
2004150901.tif 000017
【図２】
2004150901.tif 000018
【図３】
2004150901.tif 000019
【図４】
2004150901.tif 000020
【図５】
2004150901.tif 000021
【図６】
2004150901.tif 000022
【図７】
2004150901.tif 000023
【図８】
2004150901.tif 000024
【図９】
2004150901.tif 000025
【図１０】
2004150901.tif 000026
【図１１】
2004150901.tif 000027
【図１２】
2004150901.tif 000028
【図１３】
2004150901.tif 000029
【図１４】
2004150901.tif 000030
【図１５】
2004150901.tif 000031