特許 6738876 補正パラメータ作成装置、補正パラメータ作成方法および補正パラメータ作成プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6738876

(24)【登録日】2020年7月22日

(45)【発行日】2020年8月12日

(54)【発明の名称】補正パラメータ作成装置、補正パラメータ作成方法および補正パラメータ作成プログラム

(51)【国際特許分類】

   G01S 19/40 20100101AFI20200730BHJP

   G01C 15/00 20060101ALI20200730BHJP

【ＦＩ】

   G01S19/40

   G01C15/00 102C

【請求項の数】4

【全頁数】22

(21)【出願番号】特願2018-202894(P2018-202894)

(22)【出願日】2018年10月29日

(65)【公開番号】特開2020-71044(P2020-71044A)

(43)【公開日】2020年5月7日

【審査請求日】2019年2月5日

(73)【特許権者】

【識別番号】593154436

【氏名又は名称】アイサンテクノロジー株式会社

(74)【代理人】

【識別番号】110000578

【氏名又は名称】名古屋国際特許業務法人

(72)【発明者】

【氏名】中根 勝見

(72)【発明者】

【氏名】細井 幹広

(72)【発明者】

【氏名】市川 裕貴

【審査官】 田中 純

(56)【参考文献】

【文献】 特開２００８−１６４５０８（ＪＰ，Ａ）

【文献】 特開２００７−０９３５１３（ＪＰ，Ａ）

【文献】 特開２０１６−００８８７６（ＪＰ，Ａ）

【文献】 特開２０１７−０１５４５８（ＪＰ，Ａ）

【文献】 特開２０１４−０３５２１９（ＪＰ，Ａ）

【文献】 特開平１０−１００７３８（ＪＰ，Ａ）

【文献】 特開２００６−２５８４６１（ＪＰ，Ａ）

【文献】 特開平１１−１０２１５６（ＪＰ，Ａ）

【文献】 特開２００２−２４３８３３（ＪＰ，Ａ）

【文献】 特開２００６−２５９２０７（ＪＰ，Ａ）

【文献】 特開２０１０−１３９３３５（ＪＰ，Ａ）

【文献】 特開２０１１−１０２７３６（ＪＰ，Ａ）

【文献】 特開２０１６−０８０５８９（ＪＰ，Ａ）

【文献】 特開２０１８−０８１２５２（ＪＰ，Ａ）

【文献】 特開２０１９−０２０２６０（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１８／０２１６９３４（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｓ ５／００ − Ｇ０１Ｓ ５／１４

Ｇ０１Ｓ １９／００ − Ｇ０１Ｓ １９／５５

Ｇ０１Ｃ １／００ − Ｇ０１Ｃ １／１４

Ｇ０１Ｃ ５／００ − Ｇ０１Ｃ １５／１４

Ｇ０１Ｃ ２１／００ − Ｇ０１Ｃ ２１／３６

Ｇ０１Ｃ ２３／００ − Ｇ０１Ｃ ２５／００

(57)【特許請求の範囲】

【請求項1】

単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する補正パラメータ作成装置であって、
前記補正パラメータを作成する基準となる日である基準日を設定するように構成された基準日設定部と、
前記補正パラメータが有効な期間であるパラメータ有効期間を設定するように構成された有効期間設定部と、
前記パラメータ有効期間内に、前記パラメータ有効期間内の任意の日である使用日を設定するように構成された使用日設定部と、
前記補正パラメータを作成する対象となる地域である対象地域内に存在する複数の基準点のそれぞれについて、前記基準日を含む一定期間内の前記基準点の座標値に基づき、前記基準日における前記基準点の位置を基準日位置として算出するように構成された基準日位置算出部と、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、算出時点に得られる一定期間内の前記基準点の座標値に基づき、前記使用日における前記基準点の位置を使用日位置として算出するように構成された使用日位置算出部と、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、前記基準日位置算出部により算出された前記基準日位置と、前記使用日位置算出部により算出された前記使用日位置との差を、地殻変動量として算出するように構成された変動量算出部と、
前記対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、空間補間法を用いて、前記地点を囲む複数の前記基準点の前記地殻変動量に基づき、地殻変動推定量を算出し、算出した前記地殻変動推定量を前記補正パラメータとするように構成された推定量算出部と
を備える補正パラメータ作成装置。

【請求項2】

請求項１に記載の補正パラメータ作成装置であって、
前記使用日設定部は、前記パラメータ有効期間の中間に前記使用日を設定する補正パラメータ作成装置。

【請求項3】

単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する補正パラメータ作成方法であって、
前記補正パラメータを作成する基準となる日である基準日を設定する基準日設定手順と、
前記補正パラメータが有効な期間であるパラメータ有効期間を設定する有効期間設定手順と、
前記パラメータ有効期間内に、前記パラメータ有効期間内の任意の日である使用日を設定する使用日設定手順と、
前記補正パラメータを作成する対象となる地域である対象地域内に存在する複数の基準点のそれぞれについて、前記基準日を含む一定期間内の前記基準点の座標値に基づき、前記基準日における前記基準点の位置を基準日位置として算出する基準日位置算出手順と、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、算出時点に得られる一定期間内の前記基準点の座標値に基づき、前記使用日における前記基準点の位置を使用日位置として算出する使用日位置算出手順と、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、前記基準日位置算出手順により算出された前記基準日位置と、前記使用日位置算出手順により算出された前記使用日位置との差を、地殻変動量として算出する変動量算出手順と、
前記対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、空間補間法を用いて、前記地点を囲む複数の前記基準点の前記地殻変動量に基づき、地殻変動推定量を算出し、算出した前記地殻変動推定量を前記補正パラメータとする推定量算出手順と
を備える補正パラメータ作成方法。

【請求項4】

単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成するために、コンピュータを、
前記補正パラメータを作成する基準となる日である基準日を設定するように構成された基準日設定部、
前記補正パラメータが有効な期間であるパラメータ有効期間を設定するように構成された有効期間設定部、
前記パラメータ有効期間内に、前記パラメータ有効期間内の任意の日である使用日を設定するように構成された使用日設定部、
前記補正パラメータを作成する対象となる地域である対象地域内に存在する複数の基準点のそれぞれについて、前記基準日を含む一定期間内の前記基準点の座標値に基づき、前記基準日における前記基準点の位置を基準日位置として算出するように構成された基準日位置算出部、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、算出時点に得られる一定期間内の前記基準点の座標値に基づき、前記使用日における前記基準点の位置を使用日位置として算出するように構成された使用日位置算出部、
前記対象地域内に存在する複数の前記基準点のそれぞれについて、前記基準日位置算出部により算出された前記基準日位置と、前記使用日位置算出部により算出された前記使用日位置との差を、地殻変動量として算出するように構成された変動量算出部、及び、
前記対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、空間補間法を用いて、前記地点を囲む複数の前記基準点の前記地殻変動量に基づき、地殻変動推定量を算出し、算出した前記地殻変動推定量を前記補正パラメータとするように構成された推定量算出部
として機能させるための補正パラメータ作成プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本開示は、測定位置を地殻変動に基づいて補正するための補正パラメータを作成する補正パラメータ作成装置、方法およびプログラムに関する。

【背景技術】

【0002】

特許文献１および特許文献２には、全地球航法衛星システムを用いて位置情報を取得し、取得した位置情報が示す位置を地図上に表示する技術が知られている。
衛星測位分野では、大きく分けて、基準点の不要な単独測位と、基準点が必要な相対測位とが行われおり、一般的には後者の方が高精度である。実用化されている高精度な測位方法としては、相対測位であるＲＴＫ法があり、測量分野をはじめとして様々な分野で利用されている。ＲＴＫは、Real-Time Kinematicの略である。

【0003】

しかし近年では、ＰＰＰ、ＰＰＰ−ＲＴＫ、ＰＰＰ−ＡＲなどの高精度な単独測位が発展してきている。ＰＰＰは、Precise Point Positioningの略である。ＡＲは、Ambiguity Resolutionの略である。

【0004】

日本においては、２０１８年度より実運用が開始される準天頂衛星の補強信号として採用されている、センチメータ級測位補強サービス（ＣＬＡＳ）と、ＭＡＤＯＣＡとがその一例である。ＣＬＡＳは、Centimeter-Level Augmentation Serviceの略である。ＭＡＤＯＣＡは、Multi-GNSS Advanced Demonstration tool for Orbit and Clock Analysisの略である。高精度単独測位は、基準局を必要としないため、自動走行などの先端分野から期待が高まっている。なお、基準局が必要な場合には、基準局までの距離（すなわち、基線長）、および、基準局の切り替えなどの問題がある。

【0005】

地図の座標値は、測量成果に基づいており、時間情報である測量成果の元期に紐づけられている。このような座標を元期座標という。一方、単独測位によって得られる座標値は、観測日に紐づけられている。このような座標を今期座標という。一般的に、元期座標と今期座標とで、紐づけられた時間は異なる。つまり、時間経過に伴う影響、特に地殻変動の影響を大きく受ける。よって、単独測位によって得られた座標値と地図の座標値とは、時間経過とともに乖離するという問題が生じる。

【0006】

上記の問題の解決方法として、元期座標から今期座標の間に生じる地殻変動による移動量を補正する方法が考えられる。実例として、日本においては、国土地理院が公開している地殻変動補正パラメータ（以下、ＳＤパラメータ）を使用して地殻変動量を補正している。この補正方法をセミ・ダイナミック補正といい、地殻変動による位置の乖離を解消し、地図と整合性のとれた座標値を得ることができる。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】特開２００１−０５９７２６号公報

【特許文献2】特開２００４−１４４７１１号公報

【発明の概要】

【発明が解決しようとする課題】

【0008】

しかし、ＳＤパラメータの更新頻度は１年である。そして、地殻変動が活発な場所では年間約１０ｃｍも移動するため、高精度測位の精度よりも、１年間での地殻変動量の方が大きくなってしまう。また、ＳＤパラメータは毎年１月１日の位置情報を利用しているため、地殻変動補正の効果が時間経過とともに単調に減少する。
また、作成時点の基準点の位置から計算したパラメータは提供時には既に劣化し、使用時に最適に補正されない。また、基準点座標を連続観測していると、なかには地殻変動に依らないノイズ（システム雑音や観測雑音）を含む座標値が発生するので、それがパラメータに影響することがある。

【0009】

本開示は、地殻変動補正の効果を向上させることを目的とする。

【課題を解決するための手段】

【0010】

本開示の一態様は、単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する補正パラメータ作成装置である。
そして、本開示の補正パラメータ作成装置は、基準日設定部と、有効期間設定部と、使用日設定部と、基準日位置算出部と、使用日位置算出部と、変動量算出部と、推定量算出部とを備える。ここで、基準日とは、地図等の位置の基準になる日を意味する。パラメータ有効期間とは、補正パラメータが有効な期間を意味する。使用日とは、パラメータ有効期間内の任意の日を意味する。

【0011】

基準日設定部は、基準日を設定するように構成される。有効期間設定部は、補正パラメータが有効な期間であるパラメータ有効期間を設定するように構成される。使用日設定部は、パラメータ有効期間内に、使用日を設定するように構成される。

【0012】

基準日位置算出部は、補正パラメータを作成する対象となる地域である対象地域内に存在する複数の基準点のそれぞれについて、基準日を含む一定期間内の基準点の座標値に基づき、基準日における基準点の位置を基準日位置として算出するように構成される。

【0013】

使用日位置算出部は、対象地域内に存在する複数の基準点のそれぞれについて、算出時点に得られる一定期間内の基準点の座標値に基づき、使用日における基準点の位置を使用日位置として算出するように構成される。

【0014】

変動量算出部は、対象地域内に存在する複数の基準点のそれぞれについて、基準日位置算出部により算出された基準日位置と、使用日位置算出部により算出された使用日位置との差を、地殻変動量として算出するように構成される。

【0015】

推定量算出部は、対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、空間補間法を用いて、その地点を囲む複数の基準点の地殻変動量に基づき、地殻変動推定量を算出し、算出した地殻変動推定量を補正パラメータとするように構成される。

【0016】

このように構成された本開示の補正パラメータ作成装置は、パラメータ有効期間内の使用日の地殻変動量に基づいて、補正パラメータを作成する。このため、本開示の補正パラメータ作成装置は、補正パラメータによる地殻変動補正の効果が、パラメータ有効期間内で最大となるようにすることができる。すなわち、本開示の補正パラメータ作成装置は、パラメータ有効期間の始期から使用日までは、地殻変動補正の効果が単調に増加し、使用日からパラメータ有効期間の終期までは、地殻変動補正の効果が単調に減少するようにすることが期待できる。これにより、本開示の補正パラメータ作成装置は、パラメータ有効期間全体として地殻変動補正の効果を向上させることができる。

【0017】

さらに、本開示の補正パラメータ作成装置は、一定期間内の基準点の座標値に基づき、基準日および使用日における基準点の位置（すなわち、基準日位置および使用日位置）を算出して、補正パラメータを作成する。これにより、本開示の補正パラメータ作成装置は、補正パラメータに対して、地殻変動によらないノイズによる基準点位置の変動の影響が及ぶのを抑制することができる。

【0018】

本開示の一態様では、使用日設定部は、パラメータ有効期間の中間に使用日を設定するようにしてもよい。これにより、本開示の補正パラメータ作成装置は、パラメータ有効期間内で最適に補正パラメータによる補正をさせることが期待できる。

【0019】

本開示の別の態様は、単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する補正パラメータ作成方法である。

【0020】

そして、本開示の補正パラメータ作成方法は、基準日設定手順と、有効期間設定手順と、使用日設定手順と、基準日位置算出手順と、使用日位置算出手順と、変動量算出手順と、推定量算出手順とを備える。

【0021】

基準日設定手順は、基準日を設定する。有効期間設定手順は、補正パラメータが有効な期間であるパラメータ有効期間を設定する。使用日設定手順は、パラメータ有効期間内に、使用日を設定する。

【0022】

基準日位置算出手順は、補正パラメータを作成する対象となる地域である対象地域内に存在する複数の基準点のそれぞれについて、基準日を含む一定期間内の基準点の座標値に基づき、基準日における基準点の位置を基準日位置として算出する。

【0023】

使用日位置算出手順は、対象地域内に存在する複数の基準点のそれぞれについて、算出時点に得られる一定期間内の基準点の座標値に基づき、使用日における基準点の位置を使用日位置として算出する。

【0024】

変動量算出手順は、対象地域内に存在する複数の基準点のそれぞれについて、基準日位置算出手順により算出された基準日位置と、使用日位置算出手順により算出された使用日位置との差を、地殻変動量として算出する。

【0025】

推定量算出手順は、対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、空間補間法を用いて、その地点を囲む複数の基準点の地殻変動量に基づき、地殻変動推定量を算出し、算出した地殻変動推定量を補正パラメータとする。

【0026】

本開示の補正パラメータ作成方法は、本開示の補正パラメータ作成装置にて実行される方法であり、当該方法を実行することで、本開示の補正パラメータ作成装置と同様の効果を得ることができる。

【0027】

本開示の更に別の態様は、コンピュータを、基準日設定部、有効期間設定部、使用日設定部、基準日位置算出部、使用日位置算出部、変動量算出部、及び、推定量算出部として機能させるための補正パラメータ作成プログラムである。

【0028】

本開示の補正パラメータ作成プログラムによって制御されるコンピュータは、本開示の補正パラメータ作成装置の一部を構成することができ、本開示の補正パラメータ作成装置と同様の効果を得ることができる。

【図面の簡単な説明】

【0029】

【図1】補正パラメータ作成装置の構成を示すブロック図である。

【図2】第１実施形態の補正パラメータ作成処理を示すフローチャートである。

【図3】有効期間中間日の算出方法を説明する図である。

【図4】第１実施形態においてパラメータ基準日での基準点の座標を算出する方法を説明する図である。

【図5】第１実施形態において有効期間中間日での基準点の座標を算出する方法を説明する図である。

【図6】地殻変動量ベクトルを算出する方法を説明する図である。

【図7】経験バリオグラムを示す図である。

【図8】第２実施形態の補正パラメータ作成処理を示すフローチャートである。

【図9】ＳＤパラメータの線形補間方法を説明する図である。

【図10】第２実施形態においてパラメータ基準日での基準点の座標を算出する方法を説明する図である。

【図11】第２実施形態において有効期間中間日での基準点の座標を算出する方法を説明する図である。

【図12】有効期間中間日における地殻変動量を算出する方法を説明する図である。

【発明を実施するための形態】

【0030】

（第１実施形態）
以下に本開示の第１実施形態を図面とともに説明する。
本実施形態の補正パラメータ作成装置１は、図１に示すように、表示部１１と、操作入力部１２と、データ記憶部１３と、データ入出力部１４と、制御部１５とを備える。

【0031】

表示部１１は、図示しない表示装置を備え、表示装置の表示画面に各種画像を表示する。
操作入力部１２は、図示しないキーボードおよびマウスを介して使用者が行った入力操作を特定するための入力操作情報を出力する。

【0032】

データ記憶部１３は、各種データを記憶するための記憶装置である。
データ入出力部１４は、有線または無線で接続された外部機器との間でデータの入出力を行う。

【0033】

制御部１５は、ＣＰＵ、ＲＯＭおよびＲＡＭ等を備えたマイクロコンピュータを中心に構成される。マイクロコンピュータの各種機能は、ＣＰＵが非遷移的実体的記録媒体に格納されたプログラムを実行することにより実現される。この例では、ＲＯＭが、プログラムを格納した非遷移的実体的記録媒体に該当する。また、このプログラムの実行により、プログラムに対応する方法が実行される。なお、ＣＰＵが実行する機能の一部または全部を、一つあるいは複数のＩＣ等によりハードウェア的に構成してもよい。また、制御部１５を構成するマイクロコンピュータの数は１つでも複数でもよい。

【0034】

次に、制御部１５が実行する補正パラメータ作成処理の手順を説明する。補正パラメータ作成処理は、補正パラメータ作成処理を実行するために制御部１５に記憶された補正パラメータ作成プログラム２０を使用者の入力操作により起動することで実行される。なお、補正パラメータ作成プログラム２０は、補正パラメータ作成装置１に予めインストールされていてもよいし、記録媒体またはネットワークを介してインストールされるようにしてもよい。記録媒体としては、例えば光ディスク、磁気ディスクおよび半導体メモリなどが挙げられる。

【0035】

補正パラメータ作成処理が実行されると、制御部１５は、図２に示すように、まず、Ｓ１０にて、日付計算処理を実行する。日付計算処理では、制御部１５は、まず、有効期間開始日と有効期間日数を示す有効期間情報を入力するための画像を表示部１１の表示画面に表示する。そして、制御部１５は、有効期間情報が操作入力部１２から入力されるまで待機する。そして、有効期間情報が入力されると、制御部１５は、図３に示すように、式（１）により、有効期間中間日を算出する。式（１）におけるＤＭは有効期間中間日、ＤＢは有効期間開始日、ＤＵＲは有効期間日数である。また、式（１）における［＊］は、ガウス記号である。なお、ガウス記号内の＊は数字である。

【0036】

【数1】

Ｓ１０の処理が終了すると、図２に示すように、Ｓ２０にて、制御部１５は、メッシュ取込処理を実行する。メッシュ取込処理では、制御部１５は、まず、補正パラメータを作成する対象となる地域（以下、対象地域）を示す対象地域情報を入力するための画像を表示部１１の表示画面に表示する。そして制御部１５は、対象地域情報が操作入力部１２から入力されるまで待機する。そして、対象地域情報が入力されると、制御部１５は、対象地域情報のメッシュデータをデータ記憶部１３から取得する。メッシュは、対象地域を緯度および経度に基づいて格子状に区切って形成される。そしてメッシュデータは、格子状に形成されるメッシュの格子点（以下、メッシュ頂点）における緯度および経度を示すデータである。なお、本実施形態では、メッシュデータは、補正パラメータ作成処理を開始する前に予めデータ記憶部１３に記憶される。

【0037】

またＳ３０にて、制御部１５は、位置情報取得処理を実行する。位置情報取得処理では、制御部１５は、対象地域内に存在するＣＯＲＳの位置を示す基準点位置情報をデータ記憶部１３から取得する。ＣＯＲＳは、測量における基準点であり、測位衛星からの信号を受信することにより自身の位置が計測される。ＣＯＲＳは、Continuously Operating Reference Stationの略である。基準点位置情報は、少なくとも、ＣＯＲＳの緯度、経度、楕円体高、年月日および時刻を示す情報である。

【0038】

位置情報取得処理では、具体的に、制御部１５は、まず、後述するパラメータ基準日と、基準点位置情報の数（以下、位置取得数）とを示す取得情報を入力するための画像を表示部１１の表示画面に表示する。そして制御部１５は、取得情報が操作入力部１２から入力されるまで待機する。そして、取得情報が入力されると、制御部１５は、複数のＣＯＲＳ毎に、パラメータ基準日を中心とした位置取得数分の基準点位置情報を取得する。さらに制御部１５は、複数のＣＯＲＳ毎に、年月日および時刻が新しい順に位置取得数分の基準点位置情報を取得する。位置取得数は、例えば、１００〜２００である。

【0039】

次にＳ４０にて、制御部１５は、位置決定処理を実行する。位置決定処理では、制御部１５は、Ｓ３０で基準点位置情報を取得した複数のＣＯＲＳ毎に、パラメータ基準日におけるＣＯＲＳの位置と、有効期間中間日におけるＣＯＲＳの位置とを算出する。

【0040】

位置決定処理では、具体的には、制御部１５は、まず、複数のＣＯＲＳのうち、パラメータ基準日および有効期間中間日における位置が算出されていない１つのＣＯＲＳを選択する。そして制御部１５は、選択したＣＯＲＳにおいて、パラメータ基準日を中心とした位置取得数分の基準点位置情報の中から、最新の年月日および時刻の基準点位置情報（以下、第１最新基準点位置情報）と、最古の年月日および時刻の基準点位置情報（以下、第１最古基準点位置情報）とを抽出する。さらに制御部１５は、第１最古基準点位置情報における年月日および時刻を第１算出期間開始時ＤＵＰＢ１とし、第１最新基準点位置情報における年月日および時刻を第１算出期間終了時ＤＵＰＥ１とする。なお、第１算出期間開始時ＤＵＰＢ１から第１算出期間終了時ＤＵＰＥ１までの期間を第１算出期間ＤＵＰ１という。

【0041】

そして制御部１５は、第１算出期間ＤＵＰ１内の基準点位置情報を用い、緯度、経度および楕円体高のそれぞれについて線形回帰直線を算出する。さらに制御部１５は、緯度、経度および楕円体高のそれぞれで算出された線形回帰直線について、パラメータ基準日における緯度、経度および楕円体高を算出する。図４に示すように、点Ｐ１，Ｐ２，Ｐ３，Ｐ４，Ｐ５，Ｐ６，Ｐ７，Ｐ８は、第１算出期間ＤＵＰ１内の基準点位置情報における緯度、経度または楕円体高を示す。直線Ｌ１は、線形回帰直線を示す。点Ｐ９は、線形回帰直線に基づいて算出されたパラメータ基準日ＤＲにおける緯度、経度または楕円体高である。なお、線形回帰直線の代わりに、状態空間推定やディープラーニングを利用したリカレントネットワーク等の時系列推定を利用してもよい。

【0042】

次に制御部１５は、選択したＣＯＲＳにおいて、年月日および時刻が新しい順に取得した位置取得数分の基準点位置情報の中から、最新の年月日および時刻の基準点位置情報（以下、第２最新基準点位置情報）と、最古の年月日および時刻の基準点位置情報（以下、第２最古基準点位置情報）とを抽出する。さらに制御部１５は、第２最古基準点位置情報における年月日および時刻を第２算出期間開始時ＤＵＰＢ２とし、第２最新基準点位置情報における年月日および時刻を第２算出期間終了時ＤＵＰＥ２とする。なお、第２算出期間開始時ＤＵＰＢ２から第２算出期間終了時ＤＵＰＥ２までの期間を第２算出期間ＤＵＰ２という。

【0043】

そして制御部１５は、第２算出期間ＤＵＰ２内の基準点位置情報を用い、緯度、経度および楕円体高のそれぞれについて線形回帰直線を算出する。そして、緯度、経度および楕円体高のそれぞれで算出された線形回帰直線について、有効期間中間日における緯度、経度および楕円体高を算出する。図５に示すように、点Ｐ１１，Ｐ１２，Ｐ１３，Ｐ１４，Ｐ１５は、第２算出期間ＤＵＰ２内の基準点位置情報における緯度、経度または楕円体高を示す。直線Ｌ２は、線形回帰直線を示す。点Ｐ１６は、線形回帰直線に基づいて算出された有効期間中間日ＤＭにおける緯度、経度または楕円体高である。

【0044】

そして制御部１５は、パラメータ基準日および有効期間中間日における位置が算出されていないＣＯＲＳが存在しているか否かを判断する。ここで、位置が算出されていないＣＯＲＳが存在している場合には、ＣＯＲＳを新たに選択して、上述の処理を繰り返す。一方、位置が算出されていないＣＯＲＳが存在しない場合には、位置決定処理を終了する。

【0045】

Ｓ４０の処理が終了すると、図２に示すように、制御部１５は、Ｓ５０にて、変動量算出処理を実行する。具体的には、制御部１５は、図６に示すように、Ｓ３０で基準点位置情報を取得した複数のＣＯＲＳ毎に、式（２）により、パラメータ基準日におけるＣＯＲＳの座標値ＣＲ（すなわち、緯度、経度および楕円体高）を始点とし、有効期間中間日におけるＣＯＲＳの座標値ＣＭを終点とする地殻変動量ベクトルＶＭを算出する。

【0046】

【数2】

Ｓ５０の処理が終了すると、図２に示すように、制御部１５は、Ｓ６０にて、経験バリオグラム作成処理を実行する。

【0047】

経験バリオグラム作成処理では、具体的に、制御部１５は、まず、ビン幅を示すビン幅情報を入力するための画像を表示部１１の表示画面に表示する。そして制御部１５は、ビン幅情報が操作入力部１２から入力されるまで待機する。そして、ビン幅情報が入力されると、制御部１５は、Ｓ３０で基準点位置情報を取得した複数のＣＯＲＳの中から、後述する非類似度が算出されていない２つのＣＯＲＳの組み合わせ（以下、ＣＯＲＳ対）を選択する。そして制御部１５は、選択したＣＯＲＳ対において、二点間の距離（以下、点間距離）を算出する。２つのＣＯＲＳ間の距離として、平面投影した座標値や地心座標系における直線距離、測地線長を利用してもよい。

【0048】

そして制御部１５は、選択したＣＯＲＳ対における各成分（すなわち、緯度、経度および楕円体高）の非類似度ＤＩＳを、式（３）により算出する。式（３）におけるＳＭ１は、選択したＣＯＲＳ対の一方における地殻変動量ベクトルＶＭの各成分である。式（３）におけるＳＭ２は、選択したＣＯＲＳ対の他方における地殻変動量ベクトルＶＭの各成分である。

【0049】

【数3】

そして制御部１５は、図７に示すように、地殻変動量ベクトルＶＭの各成分について、点間距離を横軸とし、非類似度ＤＩＳを縦軸とした二次元直交座標系に、算出した点間距離と非類似度ＤＩＳとを示す座標点を配置することにより、経験バリオグラムを作成する。図７は、点間距離を示す横軸をビン幅で分割したＢＮ１〜ＢＮ７内に座標点Ｐ２１〜Ｐ３２が配置された経験バリオグラムを示している。

【0050】

次に制御部１５は、非類似度ＤＩＳが算出されていないＣＯＲＳ対が存在しているか否かを判断する。ここで、非類似度ＤＩＳが算出されていないＣＯＲＳ対が存在している場合には、ＣＯＲＳ対を新たに選択して、上述の処理を繰り返す。一方、非類似度ＤＩＳが算出されていないＣＯＲＳ対が存在しない場合には、経験バリオグラム作成処理を終了する。

【0051】

Ｓ６０の処理が終了すると、図２に示すように、制御部１５は、Ｓ７０にて、モデル推定処理を実行する。モデル推定処理では、制御部１５は、経験バリオグラムから、式（４）で表されるバリオグラム関数γ（ｒ）を推定する。式（４）におけるＣ_３は、式（５）で表される。

【0052】

モデル推定処理では、具体的に、制御部１５は、まず、Ｓ６０で作成された経験バリオグラムを用いて、推定用データ列（ｘ_ｉ，ｙ_ｉ）を作成する。ｘ_ｉは、ビンｉの階級の中心である。ｙ_ｉは、ビンｉに含まれるＮｉ個の非類似度ＤＩＳの平均値である。

【0053】

例えば、図７において、ビンＢＮ１，ＢＮ２，ＢＮ３，ＢＮ４，ＢＮ５，ＢＮ６，ＢＮ７はそれぞれ、ビン１，２，３，４，５，６，７に相当する。そして、ビン１に相当するビンＢＮ１内には、座標点Ｐ２１に対応する１個の非類似度ＤＩＳが含まれる。また、ビン５に相当するビンＢＮ５内には、座標点Ｐ２７，Ｐ２８に対応する２個の非類似度ＤＩＳが含まれる。

【0054】

ｘ_ｉ，ｙ_ｉはそれぞれ、式（６），（７）で表される。式（６）におけるＢＷは、ビン幅である。

【0055】

【数4】

そして制御部１５は、最小二乗法によって、作成された複数の推定用データ列（ｘ_ｉ，ｙ_ｉ）をバリオグラム関数γ（ｒ）で近似することにより、バリオグラム関数γ（ｒ）のパラメータａ，ｂを算出する。

【0056】

制御部１５は、上記のパラメータａ，ｂの算出を、地殻変動量ベクトルＶＭの各成分について実行する。
Ｓ７０の処理が終了すると、図２に示すように、制御部１５は、Ｓ８０にて、クリギング法計算処理を実行する。クリギング法計算処理では、具体的に、制御部１５は、まず、通常型クリギング法の係数行列Ａを算出する。係数行列Ａは、式（８）で表される。式（８）におけるｎは、Ｓ３０で基準点位置情報を取得したＣＯＲＳの数である。式（８）におけるｒ_ｉ，ｊは、ｉで特定されるＣＯＲＳとｊで特定されるＣＯＲＳとの点間距離である。ｉ，ｊは、１〜ｎの整数である。クリギング法は、空間補間法の一種であり、国土地理院が提供しているＳＤパラメータの算出にも使用されている。

【0057】

【数5】

そして制御部１５は、係数行列ＡのＬＵ分解を行う。すなわち、制御部１５は、係数行列Ａを、下三角行列Ｌと、上三角行列Ｕとの積に分解する。なお、ＬＵ分解のアルゴリズムとして、Ｃｒｏｕｔ法およびＤｏｏｌｉｔｔｌｅ法が知られている。

【0058】

制御部１５は、上記の係数行列Ａの算出と係数行列ＡのＬＵ分解とを、地殻変動量ベクトルＶＭの各成分について実行する。
Ｓ８０の処理が終了すると、制御部１５は、Ｓ９０にて、地殻変動量推定処理を実行する。

【0059】

地殻変動量推定処理では、具体的に、制御部１５は、まず、Ｓ２０にてメッシュデータが取得された複数のメッシュ頂点の中から、後述する地殻変動推定量が算出されていないメッシュ頂点を一つ選択する。そして制御部１５は、選択したメッシュ頂点のメッシュデータ（すなわち、メッシュ頂点における緯度および経度）を取得する。次に制御部１５は、係数ベクトルｖを算出する。係数ベクトルｖは、ｖ_１，ｖ_２，ｖ_３，・・・・，ｖ_ｎ，ｖ_ｎ＋１を成分とするベクトルである。ｎは、Ｓ３０で基準点位置情報を取得したＣＯＲＳの数である。ｉ＝１，２，３，・・・，ｎである場合において、ｖ_ｉ＝γ（ｒ_ｉ）である。ｒ_ｉは、選択したメッシュ頂点と、ｉで特定されるＣＯＲＳとの間の距離である。また、ｖ_ｎ＋１＝１である。

【0060】

そして制御部１５は、変数ベクトルｗを用いて、式（９）に示す方程式を解く。変数ベクトルｗは、式（９）に示すように、ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λを成分とするベクトルである。λは、ラグランジュの未定乗数である。すなわち、制御部１５は、ラグランジュの未定乗数法により、式（９）を解く。なお、制御部１５は、式（９）を解く際に、係数行列ＡのＬＵ分解を行う。したがって、制御部１５は、Ａｗ＝（ＬＵ）ｗ＝ｖを解くことにより、ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λを算出する。そして制御部１５は、上記のｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λの算出を各成分（すなわち、緯度、経度および楕円体高）について実行する。

【0061】

【数6】

ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎを算出した後に、制御部１５は、式（１０）により、選択したメッシュ頂点の地殻変動推定量ｚを算出する。式（１０）におけるｚ_ｉは、ｉで特定されるＣＯＲＳにおける各成分（すなわち、緯度、経度および楕円体高）の地殻変動量である。

【0062】

【数7】

さらに制御部１５は、式（１１）により、メッシュ頂点における地殻変動推定量zの予測誤差σを算出する。なお、式（１１）における「・」はベクトルの内積を表す。

【0063】

【数8】

そして制御部１５は、上記の地殻変動推定量ｚおよび予測誤差σの算出を、各成分（すなわち、緯度、経度および楕円体高）について実行する。

【0064】

各成分の地殻変動推定量ｚおよび予測誤差σの算出が終了すると、制御部１５は、地殻変動推定量ｚおよび予測誤差σが算出されていないメッシュ頂点が存在しているか否かを判断する。ここで、地殻変動推定量ｚおよび予測誤差σが算出されていないメッシュ頂点が存在している場合には、メッシュ頂点を新たに選択して、上述の処理を繰り返す。一方、地殻変動推定量ｚおよび予測誤差σが算出されていないメッシュ頂点が存在しない場合には、地殻変動量推定処理を終了する。

【0065】

Ｓ９０の処理が終了すると、制御部１５は、Ｓ１００にて、地殻変動推定量出力処理を実行する。地殻変動推定量出力処理では、具体的に、制御部１５は、Ｓ９０で算出した地殻変動推定量ｚを示す地殻変動推定量情報を、補正パラメータ作成装置１の外部へ出力する。そして、Ｓ１００の処理が終了すると、制御部１５は、Ｓ１１０にて、予測誤差出力処理を実行する。予測誤差出力処理では、具体的に、制御部１５は、Ｓ９０で算出した予測誤差σを示す予測誤差情報を、補正パラメータ作成装置１の外部へ出力する。

【0066】

そして、Ｓ１１０の処理が終了すると、制御部１５は、補正パラメータ作成処理を終了する。
このように構成された補正パラメータ作成装置１は、単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する。

【0067】

そして補正パラメータ作成装置１は、パラメータ基準日を設定する。また補正パラメータ作成装置１は、有効期間を設定する。また補正パラメータ作成装置１は、有効期間中間日を設定する。ここで、パラメータ基準日とは、地図等の位置の基準になる日を意味する。有効期間とは、補正パラメータが有効な期間を意味する。

【0068】

そして補正パラメータ作成装置１は、補正パラメータを作成する対象となる対象地域内に存在する複数のＣＯＲＳのそれぞれについて、基準日を含む一定期間内のＣＯＲＳの座標値に基づき、パラメータ基準日におけるＣＯＲＳの位置（以下、基準日位置）を算出する。

【0069】

また補正パラメータ作成装置１は、対象地域内に存在する複数のＣＯＲＳのそれぞれについて、算出時点に得られる一定期間内のＣＯＲＳの座標値に基づき、有効期間中間日におけるＣＯＲＳの位置（以下、使用日位置）を算出する。

【0070】

さらに補正パラメータ作成装置１は、対象地域内に存在する複数のＣＯＲＳのそれぞれについて、算出された基準日位置と、算出された使用日位置との差を、地殻変動量ベクトルＶＭとして算出する。

【0071】

そして補正パラメータ作成装置１は、対象地域を格子状に区切って形成された複数のメッシュのそれぞれの格子点にあたる地点について、クリギング法を用いて、その地点を囲む複数のＣＯＲＳの地殻変動量ベクトルＶＭに基づき、地殻変動推定量ｚを算出し、算出した地殻変動推定量ｚを補正パラメータとする。

【0072】

このように補正パラメータ作成装置１は、有効期間内の使用日の地殻変動量に基づいて、補正パラメータを作成する。このため、補正パラメータ作成装置１は、補正パラメータによる地殻変動補正の効果が、有効期間内で最大となるようにすることができる。すなわち、補正パラメータ作成装置１は、有効期間の始期から有効期間中間日までは、地殻変動補正の効果が単調に増加し、有効期間中間日から有効期間の終期までは、地殻変動補正の効果が単調に減少するようにすることが期待できる。これにより、補正パラメータ作成装置１は、有効期間全体として地殻変動補正の効果を向上させることができる。

【0073】

さらに補正パラメータ作成装置１は、一定期間内のＣＯＲＳの座標値に基づき、パラメータ基準日および有効期間中間日におけるＣＯＲＳの位置（すなわち、基準日位置および使用日位置）を算出して、補正パラメータを作成する。これにより、補正パラメータ作成装置１は、補正パラメータに対して、地殻変動によらないノイズによるＣＯＲＳ位置の変動の影響が及ぶのを抑制することができる。

【0074】

また補正パラメータ作成装置１は、有効期間の中間に有効期間中間日を設定する。これにより、補正パラメータ作成装置１は、有効期間内で最適に補正パラメータによる補正をさせることが期待できる。

【0075】

以上説明した実施形態において、Ｓ３０は基準日設定部および基準日設定手順としての処理に相当し、Ｓ１０は有効期間設定部および有効期間設定手順としての処理に相当し、Ｓ１０は使用日設定部および使用日設定手順としての処理に相当する。

【0076】

また、Ｓ４０は基準日位置算出部および基準日位置算出手順としての処理に相当し、Ｓ４０は使用日位置算出部および使用日位置算出手順としての処理に相当し、Ｓ５０は変動量算出部および変動量算出手順としての処理に相当し、Ｓ６０〜Ｓ９０は推定量算出部および推定量算出手順としての処理に相当する。

【0077】

また、パラメータ基準日は基準日に相当し、有効期間はパラメータ有効期間に相当し、有効期間中間日は使用日に相当し、ＣＯＲＳは基準点に相当し、地殻変動量ベクトルＶＭは地殻変動量に相当する。

【0078】

（第２実施形態）
以下に本開示の第２実施形態を図面とともに説明する。なお第２実施形態では、第１実施形態と異なる部分を説明する。共通する構成については同一の符号を付す。

【0079】

第２実施形態の補正パラメータ作成装置１は、補正パラメータ作成処理が変更された点が第１実施形態と異なる。
第２実施形態の補正パラメータ作成処理が実行されると、制御部１５は、図８に示すように、まず、Ｓ２１０にて、日付計算処理を実行する。日付計算処理では、制御部１５は、まず、有効期間開始日と有効期間日数を示す有効期間情報を入力するための画像を表示部１１の表示画面に表示する。そして、制御部１５は、有効期間情報が操作入力部１２から入力されるまで待機する。そして、有効期間情報が入力されると、制御部１５は、Ｓ１０と同様にして、式（１）により、有効期間中間日を算出する。さらに制御部１５は、有効期間中間日の年における１月１日をパラメータ基準日とする。但し、制御部１５は、有効期間中間日が１月１日から３月３１日の場合は、有効期間中間日の前年における１月１日をパラメータ基準日とする。例えば、有効期間中間日が２０１８年４月１３日である場合には、パラメータ基準日は２０１８年１月１日である。

【0080】

Ｓ２１０の処理が終了すると、制御部１５は、Ｓ２２０にて、ＳＤパラメータ取得処理を実行する。
ＳＤパラメータ取得処理では、制御部１５は、まず、補正パラメータを作成する対象となる地域（以下、対象地域）を示す対象地域情報を入力するための画像を表示部１１の表示画面に表示する。そして制御部１５は、対象地域情報が操作入力部１２から入力されるまで待機する。そして、対象地域情報が入力されると、制御部１５は、対象地域情報のＳＤパラメータをデータ記憶部１３から取得する。

【0081】

ここで、ＳＤパラメータについて説明する。セミ・ダイナミック補正では、相対測位による測位結果に対して、地殻変動による歪みの影響を補正する。ＳＤパラメータは、今期座標と対象地図の測地系の元期座標との差であり、ここで、今期座標は、１年度（４月１日から次年の３月３１日まで）を通して固定のものが使用される。

【0082】

電子基準点は、国土地理院により全国約１３００箇所に設置されている。そして国土地理院は、電子基準点の観測データに基づいて、電子基準点の位置を高精度に計測している。また国土地理院は、電子基準点の日々の座標値を国土地理院のホームページで公開している。電子基準点の日々の座標値のうち、約２週間後に提供される、最終的に解析された座標値はＦ３解と呼ばれる。

【0083】

また、ＳＤパラメータは、地域メッシュコード毎に設定されている。地域メッシュは、日本国内の地域を緯度および経度に基づいて格子状に区切って形成される。そして地域メッシュコードは、格子状に形成される地域メッシュの南西の格子点（以下、地域メッシュ頂点）における緯度および経度を示すデータである。なお、本実施形態では、ＳＤパラメータおよび地域メッシュコードは、補正パラメータ作成処理を開始する前に予めデータ記憶部１３に記憶される。

【0084】

Ｓ２２０の処理が終了すると、制御部１５は、Ｓ２３０にて、補間計算処理を実行する。補間計算処理では、制御部１５は、まず、対象地域内の地域メッシュ頂点のうち、ＳＤパラメータが設定されていない地域メッシュ頂点を選択し、選択した地域メッシュ頂点の地域メッシュコードを取得する。

【0085】

そして、図９に示すように、制御部１５は、選択した地域メッシュ頂点ＭＲ１の東西に配置されている地域メッシュ頂点ＭＲ２，ＭＲ３を選択し、選択した地域メッシュ頂点ＭＲ２，ＭＲ３に設定されているＳＤパラメータを取得する。ここで、地域メッシュ頂点ＭＲ２，ＭＲ３のＳＤパラメータを取得することができた場合には、制御部１５は、地域メッシュ頂点ＭＲ２，ＭＲ３のＳＤパラメータの線形補間によって、選択した地域メッシュ頂点ＭＲ１のＳＤパラメータを算出する。

【0086】

一方、地域メッシュ頂点ＭＲ２，ＭＲ３の少なくとも一方のＳＤパラメータを取得することができなかった場合には、制御部１５は、選択した地域メッシュ頂点ＭＲ１の南北に配置されている地域メッシュ頂点ＭＲ４，ＭＲ５を選択し、選択した地域メッシュ頂点ＭＲ４，ＭＲ５に設定されているＳＤパラメータを取得する。ここで、地域メッシュ頂点ＭＲ４，ＭＲ５のＳＤパラメータを取得することができた場合には制御部１５は、地域メッシュ頂点ＭＲ４，ＭＲ５のＳＤパラメータの線形補間によって、選択した地域メッシュ頂点ＭＲ１のＳＤパラメータを算出する。一方、地域メッシュ頂点ＭＲ４，ＭＲ５の少なくとも一方のＳＤパラメータを取得することができなかった場合には、制御部１５は、地域メッシュ頂点ＭＲ１のＳＤパラメータの算出を行わない。

【0087】

次に制御部１５は、ＳＤパラメータが設定されていない地域メッシュ頂点の中から、これまでに選択されていない地域メッシュ頂点が存在しているか否かを判断する。ここで、選択されていない地域メッシュ頂点が存在している場合には、地域メッシュ頂点を新たに選択して、上述の処理を繰り返す。一方、選択されていない地域メッシュ頂点が存在しない場合には、補間計算処理を終了する。

【0088】

Ｓ２３０の処理が終了すると、図８に示すように、制御部１５は、Ｓ２４０にて、Ｆ３解取得処理を実行する。
Ｆ３解取得処理では、制御部１５は、対象地域内に存在する電子基準点のＦ３解をデータ記憶部１３から取得する。

【0089】

Ｆ３解取得処理では、具体的に、制御部１５は、まず、Ｆ３解の数（以下、Ｆ３解取得数）を示す取得情報を入力するための画像を表示部１１の表示画面に表示する。そして制御部１５は、取得情報が操作入力部１２から入力されるまで待機する。そして、取得情報が入力されると、制御部１５は、複数の電子基準点毎に、Ｓ２１０で算出されたパラメータ基準日を中心としたＦ３解取得数分のＦ３解を取得する。さらに制御部１５は、複数の電子基準点毎に、年月日および時刻が新しい順にＦ３解取得数分のＦ３解を取得する。Ｆ３解取得数は、例えば、１００〜２００である。

【0090】

次にＳ２５０にて、制御部１５は、Ｆ３解位置決定処理を実行する。Ｆ３解位置決定処理では、制御部１５は、Ｓ２４０でＦ３解を取得した複数の電子基準点毎に、パラメータ基準日における電子基準点の位置と、有効期間中間日における電子基準点の位置とを算出する。

【0091】

Ｆ３解位置決定処理では、具体的には、制御部１５は、まず、複数の電子基準点のうち、パラメータ基準日および有効期間中間日における位置が算出されていない１つの電子基準点を選択する。そして制御部１５は、選択した電子基準点において、式（１２），（１３），（１４）により、第１算出期間開始日ＤＵＦＢ１および第１算出期間終了日ＤＵＦＥ１を算出する。式（１２）〜（１４）におけるＤＲは、パラメータ基準日である。式（１２）〜（１４）におけるＮＵＦ３は、Ｆ３解取得数である。また、式（１３）は、Ｆ３解取得数が偶数である場合に用いられ、式（１４）は、Ｆ３解取得数が奇数である場合に用いられる。なお、第１算出期間開始日ＤＵＦＢ１から第１算出期間終了日ＤＵＦＥ１までの期間を第１算出期間ＤＵＦ１という。

【0092】

【数9】

そして制御部１５は、第１算出期間ＤＵＦ１内のＦ３解を用い、緯度、経度および楕円体高のそれぞれについて線形回帰直線を算出する。さらに制御部１５は、緯度、経度および楕円体高のそれぞれで算出された線形回帰直線について、パラメータ基準日における緯度、経度および楕円体高を算出する。図１０に示すように、点Ｐ４１，Ｐ４２，Ｐ４３，Ｐ４４，Ｐ４５，Ｐ４６，Ｐ４７，Ｐ４８は、第１算出期間ＤＵＦ１内のＦ３解における緯度、経度または楕円体高を示す。直線Ｌ３は、線形回帰直線を示す。点Ｐ４９は、線形回帰直線に基づいて算出されたパラメータ基準日ＤＲにおける緯度、経度または楕円体高である。

【0093】

次に制御部１５は、選択した電子基準点において、式（１５），（１６）により、第２算出期間開始日ＤＵＦＢ２および第２算出期間終了日ＤＵＦＥ２を算出する。式（１５）におけるＤＬは、Ｓ２４０で取得した最新のＦ３解の日付である。式（１６）におけるＮＵＦ３は、Ｆ３解取得数である。なお、第２算出期間開始日ＤＵＦＢ２から第２算出期間終了日ＤＵＦＥ２までの期間を第２算出期間ＤＵＦ２という。

【0094】

【数10】

そして制御部１５は、第２算出期間ＤＵＦ２内のＦ３解を用い、緯度、経度および楕円体高のそれぞれについて線形回帰直線を算出する。さらに制御部１５は、緯度、経度および楕円体高のそれぞれで算出された線形回帰直線について、有効期間中間日における緯度、経度および楕円体高を算出する。図１１に示すように、点Ｐ５１，Ｐ５２，Ｐ５３，Ｐ５４，Ｐ５５は、第２算出期間ＤＵＦ２内のＦ３解における緯度、経度または楕円体高を示す。直線Ｌ４は、線形回帰直線を示す。点Ｐ５６は、線形回帰直線に基づいて算出された有効期間中間日ＤＭにおける緯度、経度または楕円体高である。

【0095】

そして制御部１５は、パラメータ基準日および有効期間中間日における位置が算出されていない電子基準点が存在しているか否かを判断する。ここで、位置が算出されていない電子基準点が存在している場合には、電子基準点を新たに選択して、上述の処理を繰り返す。一方、位置が算出されていない電子基準点が存在しない場合には、Ｆ３解位置決定処理を終了する。

【0096】

Ｓ２５０の処理が終了すると、図８に示すように、制御部１５は、Ｓ２６０にて、Ｆ３解変動量算出処理を実行する。具体的には、制御部１５は、Ｓ５０と同様にして、Ｓ２４０でＦ３解を取得した複数の電子基準点毎に、パラメータ基準日における電子基準点の座標値（すなわち、緯度、経度および楕円体高）を始点とし、有効期間中間日における電子基準点の座標値を終点とする地殻変動量ベクトルＶＭを算出する。

【0097】

次に制御部１５は、Ｓ２７０にて、経験バリオグラム作成処理を実行する。具体的には、制御部１５は、ＣＯＲＳの代わりに電子基準点を用いる点以外はＳ６０と同様にして、経験バリオグラムを作成する。

【0098】

そして制御部１５は、Ｓ２８０にて、Ｓ７０と同様にして、モデル推定処理を実行する。さらに制御部１５は、Ｓ２９０にて、ＣＯＲＳの代わりに電子基準点を用いる点以外はＳ８０と同様にして、クリギング法計算処理を実行する。

【0099】

次に制御部１５は、Ｓ３００にて、地殻変動量推定処理を実行する。
地殻変動量推定処理では、具体的に、制御部１５は、まず、対象地域内の地域メッシュ頂点の中から、後述する地殻変動推定量が算出されていない地域メッシュ頂点を一つ選択する。そして制御部１５は、選択した地域メッシュ頂点の地域メッシュコード（すなわち、メッシュ頂点における緯度および経度）を取得する。次に制御部１５は、係数ベクトルｖを算出する。係数ベクトルｖは、ｖ_１，ｖ_２，ｖ_３，・・・・，ｖ_ｎ，ｖ_ｎ＋１を成分とするベクトルである。ｎは、Ｓ２４０でＦ３解を取得した電子基準点の数である。ｉ＝１，２，３，・・・，ｎである場合において、ｖ_ｉ＝γ（ｒ_ｉ）である。ｒ_ｉは、選択した地域メッシュ頂点と、ｉで特定される電子基準点との間の距離である。また、ｖ_ｎ＋１＝１である。

【0100】

そして制御部１５は、変数ベクトルｗを用いて、式（９）に示す方程式を解く。変数ベクトルｗは、式（９）に示すように、ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λを成分とするベクトルである。λは、ラグランジュの未定乗数である。すなわち、制御部１５は、ラグランジュの未定乗数法により、式（９）を解く。なお、制御部１５は、式（９）を解く際に、係数行列ＡのＬＵ分解を行う。したがって、制御部１５は、Ａｗ＝（ＬＵ）ｗ＝ｖを解くことにより、ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λを算出する。そして制御部１５は、上記のｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎ，λの算出を各成分（すなわち、緯度、経度および楕円体高）について実行する。

【0101】

ｗ_１，ｗ_２，ｗ_３，・・・・，ｗ_ｎを算出した後に、制御部１５は、式（１０）により、選択した地域メッシュ頂点の地殻変動推定量ｚを算出する。式（１０）におけるｚ_ｉは、ｉで特定される電子基準点における各成分（すなわち、緯度、経度および楕円体高）の地殻変動量である。

【0102】

さらに制御部１５は、式（１１）により、選択した地域メッシュ頂点における地殻変動推定量の予測誤差σを算出する。
そして制御部１５は、上記の地殻変動推定量ｚおよび予測誤差σの算出を各成分（すなわち、緯度、経度および楕円体高）について実行する。

【0103】

各成分の地殻変動推定量ｚおよび予測誤差σの算出が終了すると、制御部１５は、地殻変動推定量ｚおよび予測誤差σが算出されていない地域メッシュ頂点が存在しているか否かを判断する。ここで、地殻変動推定量ｚおよび予測誤差σが算出されていない地域メッシュ頂点が存在している場合には、地域メッシュ頂点を新たに選択して、上述の処理を繰り返す。一方、地殻変動推定量ｚおよび予測誤差σが算出されていない地域メッシュ頂点が存在しない場合には、地殻変動量推定処理を終了する。

【0104】

Ｓ３００の処理が終了すると、制御部１５は、Ｓ３１０にて、ＳＤ／Ｒパラメータ算出処理を実行する。ＳＤ／Ｒパラメータ算出処理では、制御部１５は、対象地域内の複数の地域メッシュ毎に、地殻変動推定量とＳＤパラメータとを合成して、ＳＤ／Ｒパラメータを算出する。

【0105】

ＳＤ／Ｒパラメータ算出処理では、具体的に、制御部１５は、まず、対象地域内の地域メッシュ頂点の中から、ＳＤ／Ｒパラメータが算出されていない地域メッシュ頂点を一つ選択する。そして制御部１５は、選択した地域メッシュ頂点の地域メッシュコードを取得する。次に制御部１５は、取得した地域メッシュコードに対応するＳＤパラメータを取得する。さらに制御部１５は、取得した地域メッシュコードに対応する地殻変動推定量ｚを取得する。そして制御部１５は、式（１７）により、ＳＤ／Ｒパラメータを算出する。式（１７）におけるＶＭＳＤＲは、ＳＤ／Ｒパラメータである。式（１７）におけるＶＭＳＤは、ＳＤパラメータである。式（１７）におけるＶＭＥは、地殻変動推定量ｚ（緯度、経度、楕円体高）である。

【0106】

【数11】

図１２に示すように、有効期間中間日ＤＭにおける地殻変動量は、パラメータ基準日ＤＲにおける地殻変動量（すなわち、ＳＤパラメータ）と、地殻変動推定量とを加算することにより算出される。

【0107】

次に制御部１５は、ＳＤ／Ｒパラメータが算出されていない地域メッシュ頂点が存在しているか否かを判断する。ここで、ＳＤ／Ｒパラメータが算出されていない地域メッシュ頂点が存在している場合には、地域メッシュ頂点を新たに選択して、上述の処理を繰り返す。一方、ＳＤ／Ｒパラメータが算出されていない地域メッシュ頂点が存在しない場合には、ＳＤ／Ｒパラメータ算出処理を終了する。

【0108】

Ｓ３１０の処理が終了すると、制御部１５は、図８に示すように、Ｓ３２０にて、ＳＤ／Ｒパラメータ出力処理を実行する。ＳＤ／Ｒパラメータ出力処理では、制御部１５は、Ｓ３１０で算出したＳＤ／Ｒパラメータを示すＳＤ／Ｒパラメータ情報を、補正パラメータ作成装置１の外部へ出力する。

【0109】

Ｓ３２０の処理が終了すると、制御部１５は、Ｓ３３０にて、予測誤差出力処理を実行する。予測誤差出力処理では、制御部１５は、Ｓ３００で算出した予測誤差σを示す予測誤差情報を、補正パラメータ作成装置１の外部へ出力する。

【0110】

そして、Ｓ３３０の処理が終了すると、制御部１５は、補正パラメータ作成処理を終了する。
このように構成された補正パラメータ作成装置１は、単独測位により測定された位置を地殻変動に基づいて補正するための補正パラメータを作成する。

【0111】

そして補正パラメータ作成装置１は、パラメータ基準日を設定する。また補正パラメータ作成装置１は、有効期間を設定する。また補正パラメータ作成装置１は、有効期間中間日を設定する。

【0112】

そして補正パラメータ作成装置１は、補正パラメータを作成する対象となる対象地域内に存在する複数の電子基準点のそれぞれについて、基準日を含む一定期間内の電子基準点の座標値に基づき、パラメータ基準日における電子基準点の位置（以下、基準日位置）を算出する。

【0113】

また補正パラメータ作成装置１は、対象地域内に存在する複数の電子基準点のそれぞれについて、算出時点に得られる一定期間内の電子基準点の座標値に基づき、有効期間中間日における電子基準点の位置（以下、使用日位置）を算出する。

【0114】

さらに補正パラメータ作成装置１は、対象地域内に存在する複数の電子基準点のそれぞれについて、算出された基準日位置と、算出された使用日位置との差を、地殻変動量ベクトルＶＭとして算出する。

【0115】

そして補正パラメータ作成装置１は、対象地域を格子状に区切って形成された複数の地域メッシュのそれぞれについて、クリギング法を用いて、複数の電子基準点の地殻変動量ベクトルＶＭに基づき、地域メッシュ頂点における地殻変動推定量ｚを算出し、算出した地殻変動推定量ｚを補正パラメータとする。

【0116】

このように第２実施形態の補正パラメータ作成装置１は、有効期間内の時点での地殻変動量に基づいて、補正パラメータを作成する。このため、第２実施形態の補正パラメータ作成装置１は、第１実施形態の補正パラメータ作成装置１と同様の効果を得ることができる。

【0117】

以上説明した実施形態において、Ｓ２１０は基準日設定部および基準日設定手順としての処理に相当し、Ｓ２１０は有効期間設定部および有効期間設定手順としての処理に相当し、Ｓ２１０は使用日設定部および使用日設定手順としての処理に相当する。

【0118】

また、Ｓ２５０は基準日位置算出部および基準日位置算出手順としての処理に相当し、Ｓ２５０は使用日位置算出部および使用日位置算出手順としての処理に相当し、Ｓ２６０は変動量算出部および変動量算出手順としての処理に相当し、Ｓ２７０〜Ｓ３００は推定量算出部および推定量算出手順としての処理に相当する。また、電子基準点は基準点に相当する。

【0119】

以上、本開示の一実施形態について説明したが、本開示は上記実施形態に限定されるものではなく、種々変形して実施することができる。
例えば上記実施形態では、使用日が有効期間の中間に設定される形態を示したが、有効期間の中間に限定されるものではなく、有効期間内の時点であればよい。

【0120】

上記実施形態では、クリギング法を用いて地殻変動推定量を算出する形態を示した。しかし、クリギング法以外の空間補間法を用いて地殻変動推定量を算出するようにしてもよい。

【0121】

上記実施形態では、日付を示すパラメータ基準日および有効期間中間日を設定する形態を示した。しかし、パラメータ基準日および有効期間中間日の代わりに、日付だけではなく時刻も含んだパラメータ基準日時および有効期間中間日時を設定するようにしてもよい。

【0122】

また、上記実施形態における１つの構成要素が有する機能を複数の構成要素に分担させたり、複数の構成要素が有する機能を１つの構成要素に発揮させたりしてもよい。また、上記実施形態の構成の一部を省略してもよい。また、上記実施形態の構成の少なくとも一部を、他の上記実施形態の構成に対して付加、置換等してもよい。

【0123】

上述した補正パラメータ作成装置１の他、当該補正パラメータ作成装置１を構成要素とするシステム、当該補正パラメータ作成装置１としてコンピュータを機能させるためのプログラム、このプログラムを記録した媒体、補正パラメータ作成方法など、種々の形態で本開示を実現することもできる。

【符号の説明】

【0124】

１…補正パラメータ作成装置、１５…制御部、２０…補正パラメータ作成プログラム

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】