ホーム > 特許ランキング > 三菱電機ソフトウエア株式会社
知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024128655
(43)【公開日】2024-09-24
(54)【発明の名称】画像処理装置、画像処理方法及び画像処理プログラム
(51)【国際特許分類】
G06T 1/00 20060101AFI20240913BHJP
H04N 23/60 20230101ALI20240913BHJP
G06T 7/00 20170101ALI20240913BHJP
【ＦＩ】
G06T1/00 285
H04N23/60 500
G06T7/00 640
【審査請求】未請求
【請求項の数】9
【出願形態】ＯＬ
(21)【出願番号】P 2023037751
(22)【出願日】2023-03-10
(71)【出願人】
【識別番号】591102095
【氏名又は名称】三菱電機ソフトウエア株式会社
(74)【代理人】
【識別番号】110002491
【氏名又は名称】弁理士法人クロスボーダー特許事務所
(72)【発明者】
【氏名】進藤 嘉邦
(72)【発明者】
【氏名】西村 健志
【テーマコード（参考）】
5B057
5C122
5L096
【Ｆターム（参考）】
5B057AA13
5B057AA14
5B057CA12
5B057CA16
5B057CB13
5B057CB16
5B057CD12
5B057CD20
5B057DA07
5C122DA30
5C122EA61
5C122FH04
5C122GD11
5C122HA13
5C122HA35
5C122HA75
5C122HA88
5C122HB01
5C122HB10
5L096AA09
5L096BA08
5L096CA18
5L096DA02
5L096EA14
5L096EA26
(57)【要約】
【課題】高速かつ高精度なオルソ補正を可能にする。
【解決手段】サンプル点設定部１０１は、飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定する。パラメータ算出部１０２は、サンプル点ごとに、飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出する。ＲＰＣ算出部１０３は、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出する。
【選択図】図１
【解決手段】サンプル点設定部１０１は、飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定する。パラメータ算出部１０２は、サンプル点ごとに、飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出する。ＲＰＣ算出部１０３は、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出する。
【選択図】図１
【特許請求の範囲】
【請求項1】
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定部と、
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する画像処理装置。
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する画像処理装置。
【請求項2】
前記パラメータ算出部は、
前記ゼロドップラが成立すると推定される時刻を成立推定時刻として設定し、前記成立推定時刻における前記飛行体の位置及び速度を算出し、前記成立推定時刻における前記進行方向ベクトルと前記視点ベクトルとのなす角であるドップラ角を算出し、前記ドップラ角に基づき、前記進行方向ベクトルと前記視点ベクトルとが直交状態にある場合の前記飛行体の位置と前記成立推定時刻における前記飛行体の位置との誤差を算出し、算出した前記誤差を閾値と比較し、
前記誤差が前記閾値未満であれば、前記成立推定時刻と、前記成立推定時刻における前記飛行体の位置及び速度とを、前記ゼロドップラパラメータとして採用し、
前記誤差が前記閾値以上であれば、前記成立推定時刻を補正し、補正後の成立推定時刻について前記飛行体の位置及び速度を算出し、前記ドップラ角を算出し、前記誤差を算出し、前記誤差を前記閾値と比較する動作を、前記誤差が前記閾値未満になるまで繰り返す請求項１に記載の画像処理装置。
前記ゼロドップラが成立すると推定される時刻を成立推定時刻として設定し、前記成立推定時刻における前記飛行体の位置及び速度を算出し、前記成立推定時刻における前記進行方向ベクトルと前記視点ベクトルとのなす角であるドップラ角を算出し、前記ドップラ角に基づき、前記進行方向ベクトルと前記視点ベクトルとが直交状態にある場合の前記飛行体の位置と前記成立推定時刻における前記飛行体の位置との誤差を算出し、算出した前記誤差を閾値と比較し、
前記誤差が前記閾値未満であれば、前記成立推定時刻と、前記成立推定時刻における前記飛行体の位置及び速度とを、前記ゼロドップラパラメータとして採用し、
前記誤差が前記閾値以上であれば、前記成立推定時刻を補正し、補正後の成立推定時刻について前記飛行体の位置及び速度を算出し、前記ドップラ角を算出し、前記誤差を算出し、前記誤差を前記閾値と比較する動作を、前記誤差が前記閾値未満になるまで繰り返す請求項１に記載の画像処理装置。
【請求項3】
前記パラメータ算出部は、
前記閾値としてミリメートル単位の閾値を用いる請求項２に記載の画像処理装置。
前記閾値としてミリメートル単位の閾値を用いる請求項２に記載の画像処理装置。
【請求項4】
前記パラメータ算出部は、
サンプル点ごとに、前記飛行体とサンプル点との距離であるスラントレンジを算出し、
サンプル点ごとに、前記スラントレンジと前記ゼロドップラパラメータとを用いて、ＥＣＲ（Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｒｏｔａｔｉｎｇ）座標値を算出し、
サンプル点ごとに、ＥＣＲ座標値を緯度、経度及び楕円体高へ変換し、
前記ＲＰＣ算出部は、
サンプル点ごとに、変換により得られた緯度、経度及び楕円体高を用いて、ＲＰＣを算出する請求項１に記載の画像処理装置。
サンプル点ごとに、前記飛行体とサンプル点との距離であるスラントレンジを算出し、
サンプル点ごとに、前記スラントレンジと前記ゼロドップラパラメータとを用いて、ＥＣＲ（Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｒｏｔａｔｉｎｇ）座標値を算出し、
サンプル点ごとに、ＥＣＲ座標値を緯度、経度及び楕円体高へ変換し、
前記ＲＰＣ算出部は、
サンプル点ごとに、変換により得られた緯度、経度及び楕円体高を用いて、ＲＰＣを算出する請求項１に記載の画像処理装置。
【請求項5】
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像にオルソ補正を行うオルソ補正部を有する請求項１に記載の画像処理装置。
ＲＰＣを用いて、前記撮像画像にオルソ補正を行うオルソ補正部を有する請求項１に記載の画像処理装置。
【請求項6】
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像を加工して得られた加工画像にオルソ補正を行うオルソ補正部を有する請求項１に記載の画像処理装置。
ＲＰＣを用いて、前記撮像画像を加工して得られた加工画像にオルソ補正を行うオルソ補正部を有する請求項１に記載の画像処理装置。
【請求項7】
前記サンプル点設定部は、
前記撮像画像である合成開口レーダ画像に前記複数のサンプル点を設定する請求項１に記載の画像処理装置。
前記撮像画像である合成開口レーダ画像に前記複数のサンプル点を設定する請求項１に記載の画像処理装置。
【請求項8】
コンピュータが、飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
前記コンピュータが、サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
前記コンピュータが、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とを有する画像処理方法。
前記コンピュータが、サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
前記コンピュータが、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とを有する画像処理方法。
【請求項9】
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とをコンピュータに実行させる画像処理プログラム。
サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とをコンピュータに実行させる画像処理プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、オルソ補正が行われていない撮像画像の処理に関する。
【背景技術】
【0002】
上空からの撮像画像では、フォアショートニングのために倒れ込みが生じる。従って、地形図や道路図等のいわゆる地図に撮像画像をそのまま重ね合わせても位置関係が一致しない。オルソ補正では、撮像画像を、地図と同じく、真上から見たような傾きのない画像に変換する。
オルソ補正に関して、特許文献１に開示の技術がある。
オルソ補正に関して、特許文献１に開示の技術がある。
【0003】
特許文献１の技術では、以下の処理によって合成開口レーダ画像（以下、ＳＡＲ画像という）のオルソ補正が行われる。
（１）ＤＥＭ（Ｄｉｇｉｔａｌ Ｅｌｅｖａｔｉｏｎ Ｍｏｄｅｌ）と軌道データを用いて、画素数及び向きがＳＡＲ画像と一致する「模擬ＳＡＲ画像」と「模擬ＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」を作成する。
（２）実際のＳＡＲ画像と模擬ＳＡＲ画像とのずれを画像間のマッチング処理によって計測する。
（３）上記（２）の計測結果を使用すれば、「模擬ＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」を、「実際のＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」へと変換できる。すなわち、実際のＳＡＲ画像の各画素の正確な３次元位置（緯度、経度及び楕円体高）が分かる。
（４）上記（３）を用いてオルソ補正を行う。
（１）ＤＥＭ（Ｄｉｇｉｔａｌ Ｅｌｅｖａｔｉｏｎ Ｍｏｄｅｌ）と軌道データを用いて、画素数及び向きがＳＡＲ画像と一致する「模擬ＳＡＲ画像」と「模擬ＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」を作成する。
（２）実際のＳＡＲ画像と模擬ＳＡＲ画像とのずれを画像間のマッチング処理によって計測する。
（３）上記（２）の計測結果を使用すれば、「模擬ＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」を、「実際のＳＡＲ画像の各画素の緯度、経度及び楕円体高データ」へと変換できる。すなわち、実際のＳＡＲ画像の各画素の正確な３次元位置（緯度、経度及び楕円体高）が分かる。
（４）上記（３）を用いてオルソ補正を行う。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開２００９－１８７４４０号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
上記（２）の「ずれの計測（マッチング処理）」は、上記（４）によって得られるオルソ補正画像の位置精度を決定する最も重要な要素である。このため、上記（２）の「ずれの計測（マッチング処理）」は、可能な限り精密に実施する必要がある。画像間のマッチングは、精密に実施しようとすればするほど処理時間がかかる。また、近年、ＳＡＲ画像の画素数は増加傾向であるが、画素数が多くなるほどマッチング処理に時間がかかる。
また、マッチングをいかに正確に実施しようとしても、模擬ＳＡＲ画像の源泉となるＤＥＭのメッシュの細かさよりもさらに細かい精度でずれの計測することは困難である。入手可能なＤＥＭの精度はせいぜい１ｍ程度なので、それ以上に位置精度の良いオルソ補正画像を作成することは特許文献１の技術ではできない。
このように、特許文献１の技術には「処理時間」と「精度」の課題がある。
また、マッチングをいかに正確に実施しようとしても、模擬ＳＡＲ画像の源泉となるＤＥＭのメッシュの細かさよりもさらに細かい精度でずれの計測することは困難である。入手可能なＤＥＭの精度はせいぜい１ｍ程度なので、それ以上に位置精度の良いオルソ補正画像を作成することは特許文献１の技術ではできない。
このように、特許文献１の技術には「処理時間」と「精度」の課題がある。
【0006】
本開示は、上記のような課題を解決することを主な目的とする。具体的には、本開示は、高速かつ高精度なオルソ補正を可能にすることを主な目的とする。
【課題を解決するための手段】
【0007】
本開示に係る画像処理装置は、
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定部と、
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する。
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定部と、
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する。
【発明の効果】
【0008】
本開示によれば、高速かつ高精度なオルソ補正が可能になる。
【図面の簡単な説明】
【0009】
【図1】実施の形態１に係る画像処理装置の機能構成例を示す図。
【図2】実施の形態１に係る画像処理装置のハードウェア構成例を示す図。
【図3】実施の形態１に係る画像処理装置の動作例を示すフローチャート。
【図4】実施の形態１に係るサンプル点の配置例を示す図。
【図5】実施の形態１に係るゼロドップラ位置を示す図。
【図6】実施の形態１に係るゼロドップラ位置との誤差を示す図。
【図7】実施の形態１に係る平面内の各ベクトルの関係を示す図。
【図8】実施の形態１に係る高さ方向のサンプル点の配置例を示す図。
【図9】実施の形態１に係るスラントレンジ、衛星位置及び衛星速度によるＥＣＲ座標値の算出を示す図。
【図10】実施の形態１に係る衛星位置、オフナディア角及びスラントレンジを示す図。
【図11】実施の形態２に係る画像処理装置の機能構成例を示す図。
【発明を実施するための形態】
【0010】
以下、実施の形態を図を用いて説明する。以下の実施の形態の説明及び図面において、同一の符号を付したものは、同一の部分又は相当する部分を示す。
【0011】
【0012】
本実施の形態に係る画像処理装置１００は、コンピュータである。
画像処理装置１００の動作手順は、画像処理方法に相当する。また、画像処理装置１００の動作を実現するプログラムは、画像処理プログラムに相当する。
画像処理装置１００の動作手順は、画像処理方法に相当する。また、画像処理装置１００の動作を実現するプログラムは、画像処理プログラムに相当する。
【0013】
図２に示すように、画像処理装置１００は、ハードウェアとして、プロセッサ９０１、主記憶装置９０２、補助記憶装置９０３及び通信装置９０４を備える。
また、画像処理装置１００は、図１に示すように、機能構成として、サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４を備える。サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能は、例えば、プログラムにより実現される。
補助記憶装置９０３には、サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能を実現するプログラムが記憶されている。
これらプログラムは、補助記憶装置９０３から主記憶装置９０２にロードされる。そして、プロセッサ９０１がこれらプログラムを実行して、後述するサンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の動作を行う。
図２は、プロセッサ９０１がサンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能を実現するプログラムを実行している状態を模式的に表している。
また、画像処理装置１００は、図１に示すように、機能構成として、サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４を備える。サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能は、例えば、プログラムにより実現される。
補助記憶装置９０３には、サンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能を実現するプログラムが記憶されている。
これらプログラムは、補助記憶装置９０３から主記憶装置９０２にロードされる。そして、プロセッサ９０１がこれらプログラムを実行して、後述するサンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の動作を行う。
図２は、プロセッサ９０１がサンプル点設定部１０１、パラメータ算出部１０２、ＲＰＣ算出部１０３及びオルソ補正部１０４の機能を実現するプログラムを実行している状態を模式的に表している。
【0014】
図１において、サンプル点設定部１０１は、撮像画像１５０を取得する。
撮像画像１５０は、飛行体により上空から撮像された、オルソ補正が行われていない画像である。撮像画像１５０は、例えば、飛行体である人工衛星により撮像されたＳＡＲ画像である。以下では、撮像画像１５０が人工衛星により撮像されたＳＡＲ画像であるものとする。
撮像画像１５０では、フォアショートニングのために倒れ込みが生じる。従って、地図に撮像画像１５０をそのまま重ね合わせても位置関係が一致しない。このため、撮像画像１５０に対するオルソ補正が必要である。
サンプル点設定部１０１は、撮像画像１５０に複数の仮想点である複数のサンプル点を設定する。サンプル点の詳細は後述する。
サンプル点設定部１０１により行われる処理は、サンプル点設定処理に相当する。
撮像画像１５０は、飛行体により上空から撮像された、オルソ補正が行われていない画像である。撮像画像１５０は、例えば、飛行体である人工衛星により撮像されたＳＡＲ画像である。以下では、撮像画像１５０が人工衛星により撮像されたＳＡＲ画像であるものとする。
撮像画像１５０では、フォアショートニングのために倒れ込みが生じる。従って、地図に撮像画像１５０をそのまま重ね合わせても位置関係が一致しない。このため、撮像画像１５０に対するオルソ補正が必要である。
サンプル点設定部１０１は、撮像画像１５０に複数の仮想点である複数のサンプル点を設定する。サンプル点の詳細は後述する。
サンプル点設定部１０１により行われる処理は、サンプル点設定処理に相当する。
【0015】
パラメータ算出部１０２は、サンプル点ごとに、ゼロドップラパラメータを算出する。ゼロドップラパラメータは、ゼロドップラが成立する時刻と、当該時刻での人工衛星の位置及び速度である。ゼロドップラとは、人工衛星の進行方向ベクトルと人工衛星からサンプル点への視点ベクトルとが直交する状態である。ゼロドップラの詳細は後述する。
パラメータ算出部１０２により行われる処理は、パラメータ算出処理に相当する。
パラメータ算出部１０２により行われる処理は、パラメータ算出処理に相当する。
【0016】
ＲＰＣ算出部１０３は、サンプル点ごとのゼロドップラパラメータを用いて、サンプル点ごとに緯度、経度及び楕円体高を算出し、それら緯度、経度及び楕円体高を用いてＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出する。ＲＰＣは、緯度、経度及び楕円体高からなる３次元空間から、もとの画像のピクセル番号及びライン番号からなる２次元空間への変換を表す有理多項式（分母及び分子がそれぞれ多項式）の係数である。
ＲＰＣの詳細は後述する。
ＲＰＣ算出部１０３により行われる処理は、ＲＰＣ算出処理に相当する。
ＲＰＣの詳細は後述する。
ＲＰＣ算出部１０３により行われる処理は、ＲＰＣ算出処理に相当する。
【0017】
オルソ補正部１０４は、ＲＰＣとＤＥＭ１６０を用いて、撮像画像１５０に対してオルソ補正を行う。ＲＰＣにより、３次元空間内の各点と、各点に対応する撮像画像１５０上のピクセル番号及びライン番号が分かるので、オルソ補正部１０４は、これらを用いてリサンプリング処理によってオルソ補正を行う。
そして、オルソ補正部１０４は、オルソ補正済みの撮像画像１５０である補正済画像１７０を出力する。
そして、オルソ補正部１０４は、オルソ補正済みの撮像画像１５０である補正済画像１７０を出力する。
【0018】
【0019】
ステップＳ１０１において、サンプル点設定部１０１が撮像画像１５０を取得し、また、ＲＰＣ算出用パラメータを取得する。
【0020】
次に、ステップＳ１０２において、サンプル点設定部１０１が撮像画像１５０にサンプル点を設定する。
【0021】
ステップＳ１０３～Ｓ１０９は、撮像画像１５０のピクセル番号及びライン番号ごと（つまり、サンプル点ごと）に繰り返し行われる。
【0022】
ステップＳ１０３では、パラメータ算出部１０２が、各サンプル点の緯度及び経度を算出する。
【0023】
また、ステップＳ１０４では、パラメータ算出部１０２が、緯度及び経度からＥＣＲ（Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｒｏｔａｔｉｎｇ）座標値への変換を行う。
【0024】
ステップＳ１０５では、パラメータ算出部１０２が、サンプル点ごとに、ゼロドップラでの時刻、衛星位置及び衛星速度を算出する。
ステップＳ１０５で算出される時刻、衛星位置及び衛星速度がゼロドップラパラメータである。
ステップＳ１０５で算出される時刻、衛星位置及び衛星速度がゼロドップラパラメータである。
【0025】
ステップＳ１０６では、パラメータ算出部１０２が、サンプル点ごとに、スラントレンジを算出する。スラントレンジは、人工衛星とサンプル点との距離である。
【0026】
ステップＳ１０７では、サンプル点設定部１０１が、高さ方向のサンプル点を設定する。
【0027】
ステップＳ１０８とステップＳ１０９は高さ方向で繰り返し行われる。
【0028】
ステップＳ１０８では、パラメータ算出部１０２が、スラントレンジ及びゼロドップラパラメータを用いて、サンプル点ごとにＥＣＲ座標値を算出する。
【0029】
ステップＳ１０９では、パラメータ算出部１０２が、サンプル点ごとに、ＥＣＲ座標値から緯度、経度及び楕円体高への変換を行う。
【0030】
最後に、ステップＳ１１０において、ＲＰＣ算出部１０３が、最小二乗法によりＲＰＣを算出する。つまり、ＲＰＣ算出部１０３は、ステップＳ１０９における変換により得られたサンプル点ごとの緯度、経度及び楕円体高を用いて、最小二乗法によりＲＰＣを算出する。
ＲＰＣ算出部１０３は、算出したＲＰＣをオルソ補正部１０４に出力する。
ＲＰＣ算出部１０３は、算出したＲＰＣをオルソ補正部１０４に出力する。
【0031】
次に、図３に示す各ステップでの処理の詳細を説明する。
【0032】
ステップＳ１０１では、サンプル点設定部１０１は、以下のパラメータをＲＰＣ算出用パラメータとして、撮像画像１５０から取得する。
ピクセル数
ライン数
ビームの向き
軌道（衛星位置、速度）
軌道の時刻
ピクセル番号及びライン番号から緯度及び経度への変換係数
ピクセル数
ライン数
ビームの向き
軌道（衛星位置、速度）
軌道の時刻
ピクセル番号及びライン番号から緯度及び経度への変換係数
【0033】
サンプル点設定部１０１は、「ピクセル数」及び「ライン数」は、撮像画像ファイル、例えば、ＳＡＲイメージファイルから取得する。
また、「ビームの向き」、「軌道（衛星位置、速度）」、「軌道の時刻」及び「ピクセル番号及びライン番号から緯度及び経度への変換係数」は、例えば、撮像画像ファイルに付属するメタファイルやＳＡＲリーダファイルなどから取得する。
なお、近年の人工衛星の軌道データはＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）によって得られる。また、地上処理よって高精度化したものをＳＡＲ画像に格納することがある。このため、実力的に１ｍ未満の精度を有する軌道データが撮像画像１５０に格納されている。つまり、上述の「軌道（衛星位置、速度）」として、１ｍ未満の精度の衛星位置及び速度を取得することができる。
また、「ビームの向き」、「軌道（衛星位置、速度）」、「軌道の時刻」及び「ピクセル番号及びライン番号から緯度及び経度への変換係数」は、例えば、撮像画像ファイルに付属するメタファイルやＳＡＲリーダファイルなどから取得する。
なお、近年の人工衛星の軌道データはＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）によって得られる。また、地上処理よって高精度化したものをＳＡＲ画像に格納することがある。このため、実力的に１ｍ未満の精度を有する軌道データが撮像画像１５０に格納されている。つまり、上述の「軌道（衛星位置、速度）」として、１ｍ未満の精度の衛星位置及び速度を取得することができる。
【0034】
ステップＳ１０２では、サンプル点設定部１０１は、ピクセル方向にＮｓａｍｐ＿ｐｘ個、ライン方向にＮｓａｍｐ＿ｌｎ個のサンプル点を配置する。Ｎｓａｍｐ＿ｐｘ、Ｎｓａｍｐ＿ｌｎは、ともに１０～１００程度の値とする。この場合は、空間方向のサンプル点の総数は１００個～１０，０００個程度となる。
図４は、サンプル点の配置例を示す。図４において、黒い丸がサンプル点である。
図４は、サンプル点の配置例を示す。図４において、黒い丸がサンプル点である。
【0035】
サンプル点設定部１０１は、ステップＳ１０１で取得したＲＰＣ算出用パラメータとサンプル点が配置された撮像画像１５０をパラメータ算出部１０２に出力する。
【0036】
ステップＳ１０３では、パラメータ算出部１０２が、ステップＳ１０２で設定されたピクセル番号及びライン番号を緯度及び経度に変換する。具体的には、パラメータ算出部１０２は、ＲＰＣ算出用パラメータに含まれる変換係数（「ピクセル番号及びライン番号から緯度及び経度への変換係数」）を用いて変換を行う。
なお、撮像画像１５０によっては撮像画像１５０に変換係数が格納されていない場合がある。このような場合は、パラメータ算出部１０２は、ピクセル番号及びライン番号と緯度及び経度を関連付ける情報、例えば、撮像画像１５０上の複数個所における「ピクセル番号及びライン番号」と「緯度及び経度」の組を用いて、最小二乗法によって変換係数を算出する。そして、パラメータ算出部１０２は、算出した変換係数を用いて変換を行う。
なお、撮像画像１５０によっては撮像画像１５０に変換係数が格納されていない場合がある。このような場合は、パラメータ算出部１０２は、ピクセル番号及びライン番号と緯度及び経度を関連付ける情報、例えば、撮像画像１５０上の複数個所における「ピクセル番号及びライン番号」と「緯度及び経度」の組を用いて、最小二乗法によって変換係数を算出する。そして、パラメータ算出部１０２は、算出した変換係数を用いて変換を行う。
【0037】
ステップＳ１０４では、パラメータ算出部１０２が、ステップＳ１０３で得られたサンプル点に対する緯度及び経度を、それぞれＥＣＲ座標系上の３次元座標値（Ｘ，Ｙ，Ｚ）へ変換する。
具体的には、パラメータ算出部１０２は、以下の式（１）を用いて、変換を行う。
式（１）は、参考文献１のＰ１５～Ｐ１６に記載されている。
参考文献１：飛田幹男，“世界測地系と座標変換”（日本測量協会， ２００２）
なお、式（１）において、ａは地球長半径である。ｆは扁平率である。ｅは離心率である。Ｎは卯酉線曲率半径である。Ｂは緯度である。Ｌは経度である。Ｈｅは楕円体高である。ここでは楕円体高Ｈｅはゼロとする。また、地球長半径ａ及び扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
具体的には、パラメータ算出部１０２は、以下の式（１）を用いて、変換を行う。
式（１）は、参考文献１のＰ１５～Ｐ１６に記載されている。
参考文献１：飛田幹男，“世界測地系と座標変換”（日本測量協会， ２００２）
なお、式（１）において、ａは地球長半径である。ｆは扁平率である。ｅは離心率である。Ｎは卯酉線曲率半径である。Ｂは緯度である。Ｌは経度である。Ｈｅは楕円体高である。ここでは楕円体高Ｈｅはゼロとする。また、地球長半径ａ及び扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
【0038】
【数1】
【0039】
ステップＳ１０５では、パラメータ算出部１０２が、観測点を選択する。観測点は、ゼロドップラパラメータ算出の対象となるサンプル点である。つまり、パラメータ算出部１０２は、ゼロドップラパラメータ算出の対象となるサンプル点を順に観測点として選択する。
また、パラメータ算出部１０２は、ステップＳ１０４で求めた観測点の３次元座標値（Ｘ，Ｙ，Ｚ）に対してゼロドップラパラメータを求める。前述したように、ゼロドップラパラメータは、ゼロドップラが成立する時刻と、当該時刻での人工衛星の位置及び速度である。また、ゼロドップラは、図５に示すように、人工衛星５００の進行方向ベクトルと人工衛星５００から観測点への視点ベクトルとが直交する状態である。
パラメータ算出部１０２は、ゼロドップラが成立する瞬間（図５の状態）の観測時刻、衛星位置ベクトルと衛星速度ベクトルを求める。
以下では、衛星位置ベクトルをＰの上に矢印（→）が付されている表記とＰの横に矢印（→）が付されている表記のいずれかにより示す。また、衛星速度ベクトルをＶの上に矢印（→）が付されている表記とＶの横に矢印（→）が付されている表記のいずれかにより示す。また、他のベクトルについても同様である。つまり、アルファベットの上に矢印（→）が付されている表記とアルファベットの横に矢印（→）が付されている表記は同じベクトルを意味する。
また、パラメータ算出部１０２は、ステップＳ１０４で求めた観測点の３次元座標値（Ｘ，Ｙ，Ｚ）に対してゼロドップラパラメータを求める。前述したように、ゼロドップラパラメータは、ゼロドップラが成立する時刻と、当該時刻での人工衛星の位置及び速度である。また、ゼロドップラは、図５に示すように、人工衛星５００の進行方向ベクトルと人工衛星５００から観測点への視点ベクトルとが直交する状態である。
パラメータ算出部１０２は、ゼロドップラが成立する瞬間（図５の状態）の観測時刻、衛星位置ベクトルと衛星速度ベクトルを求める。
以下では、衛星位置ベクトルをＰの上に矢印（→）が付されている表記とＰの横に矢印（→）が付されている表記のいずれかにより示す。また、衛星速度ベクトルをＶの上に矢印（→）が付されている表記とＶの横に矢印（→）が付されている表記のいずれかにより示す。また、他のベクトルについても同様である。つまり、アルファベットの上に矢印（→）が付されている表記とアルファベットの横に矢印（→）が付されている表記は同じベクトルを意味する。
【0040】
パラメータ算出部１０２は、ゼロドップラが成立すると推定される時刻を観測時刻ｔ（成立推定時刻の例）として設定する。
そして、パラメータ算出部１０２は、観測時刻ｔにおける人工衛星５００の位置及び速度を算出する。
更に、パラメータ算出部１０２は、観測時刻ｔにおける進行方向ベクトル（衛星速度ベクトルＶ→）と視点ベクトルとのなす角であるドップラ角を算出する。
また、パラメータ算出部１０２は、ドップラ角に基づき、進行方向ベクトルと視点ベクトルとが直交状態にある場合の人工衛星５００の位置（ゼロドップラ位置）と観測時刻ｔにおける人工衛星５００の位置との誤差を算出する。そして、パラメータ算出部１０２は、算出した誤差を閾値と比較する。本実施の形態では、閾値はミリメートル単位である。
誤差が閾値未満であれば、パラメータ算出部１０２は、観測時刻ｔと、観測時刻ｔにおける人工衛星５００の位置及び速度とを、ゼロドップラパラメータとして採用する。
一方、誤差が閾値以上であれば、パラメータ算出部１０２は、観測時刻ｔを補正する。そして、パラメータ算出部１０２は、補正後の観測時刻ｔについて人工衛星５００の位置及び速度を算出し、ドップラ角を算出し、誤差を算出し、誤差を閾値と比較する動作を、誤差が閾値未満になるまで繰り返す。
以上の手順の詳細を以下にて説明する。
なお、以下の説明では、人工衛星５００の軌道はほぼ円運動（少なくとも短時間では成立）として式を簡略化している。
そして、パラメータ算出部１０２は、観測時刻ｔにおける人工衛星５００の位置及び速度を算出する。
更に、パラメータ算出部１０２は、観測時刻ｔにおける進行方向ベクトル（衛星速度ベクトルＶ→）と視点ベクトルとのなす角であるドップラ角を算出する。
また、パラメータ算出部１０２は、ドップラ角に基づき、進行方向ベクトルと視点ベクトルとが直交状態にある場合の人工衛星５００の位置（ゼロドップラ位置）と観測時刻ｔにおける人工衛星５００の位置との誤差を算出する。そして、パラメータ算出部１０２は、算出した誤差を閾値と比較する。本実施の形態では、閾値はミリメートル単位である。
誤差が閾値未満であれば、パラメータ算出部１０２は、観測時刻ｔと、観測時刻ｔにおける人工衛星５００の位置及び速度とを、ゼロドップラパラメータとして採用する。
一方、誤差が閾値以上であれば、パラメータ算出部１０２は、観測時刻ｔを補正する。そして、パラメータ算出部１０２は、補正後の観測時刻ｔについて人工衛星５００の位置及び速度を算出し、ドップラ角を算出し、誤差を算出し、誤差を閾値と比較する動作を、誤差が閾値未満になるまで繰り返す。
以上の手順の詳細を以下にて説明する。
なお、以下の説明では、人工衛星５００の軌道はほぼ円運動（少なくとも短時間では成立）として式を簡略化している。
【0041】
ステップＳ１０５－１
パラメータ算出部１０２は、観測時刻ｔの初期設定を行う。つまり、パラメータ算出部１０２は、以下の式（２）のように、初期観測時刻ｔ０を設定する。
ｔ＝ｔ０ 式（２）
ここで、初期観測時刻ｔ０は真にゼロドップラが成立する瞬間である真値時刻に対して、１／４周回以内（約２０～２５分）程度の誤差があっても許容可能である。初期観測時刻ｔ０が真値時刻に近い方が早く収束する（早く真値時刻が得られる）。
パラメータ算出部１０２は、例えば撮像画像１５０に格納されている軌道データ列の中心時刻を初期観測時刻ｔ０として使用する（格納されている軌道データが、例えば、陸域観測技術衛星２号「だいち２号」（以下、ＡＬＯＳ－２）のように２８点であれば、１４点目の時刻）。
パラメータ算出部１０２は、観測時刻ｔの初期設定を行う。つまり、パラメータ算出部１０２は、以下の式（２）のように、初期観測時刻ｔ０を設定する。
ｔ＝ｔ０ 式（２）
ここで、初期観測時刻ｔ０は真にゼロドップラが成立する瞬間である真値時刻に対して、１／４周回以内（約２０～２５分）程度の誤差があっても許容可能である。初期観測時刻ｔ０が真値時刻に近い方が早く収束する（早く真値時刻が得られる）。
パラメータ算出部１０２は、例えば撮像画像１５０に格納されている軌道データ列の中心時刻を初期観測時刻ｔ０として使用する（格納されている軌道データが、例えば、陸域観測技術衛星２号「だいち２号」（以下、ＡＬＯＳ－２）のように２８点であれば、１４点目の時刻）。
【0042】
ステップＳ１０５－２
パラメータ算出部１０２は、観測時刻ｔ（１回目であれば初期観測時刻ｔ０）における衛星位置ベクトルＰ→と衛星速度ベクトルＶ→を算出する。
具体的には、パラメータ算出部１０２は、観測時刻ｔにおける衛星位置ベクトルＰ→、衛星速度ベクトルＶ→を補間処理によって求める。
パラメータ算出部１０２は、補間方法として、例えば最小二乗法による多項式を使用することができる。撮像画像１５０に格納されている軌道データは、例えばＡＬＯＳ－２の場合は１分間隔で２８点格納されている。パラメータ算出部１０２は、これを用いて、衛星位置ベクトルＰ→の３成分（ｘ，ｙ，ｚ）それぞれを時刻を変数とした２次ないし２次以上のＮ次多項式で表現する（多項式の係数決定に最小二乗法を使う）。そして、それら３つのＮ次多項式に観測時刻ｔを与える。このようにすることで、観測時刻ｔにおける衛星位置ベクトルＰ→が求まる。衛星速度ベクトルＶ→については上記３つのＮ次多項式をそれぞれ時刻で微分した式（３つのＮ－１次多項式）を用意しておく。そして、それらＮ－１次多項式に観測時刻ｔを与える。このようにすることで、観測時刻ｔにおける衛星速度ベクトルＶ→が求まる。
パラメータ算出部１０２は、観測時刻ｔ（１回目であれば初期観測時刻ｔ０）における衛星位置ベクトルＰ→と衛星速度ベクトルＶ→を算出する。
具体的には、パラメータ算出部１０２は、観測時刻ｔにおける衛星位置ベクトルＰ→、衛星速度ベクトルＶ→を補間処理によって求める。
パラメータ算出部１０２は、補間方法として、例えば最小二乗法による多項式を使用することができる。撮像画像１５０に格納されている軌道データは、例えばＡＬＯＳ－２の場合は１分間隔で２８点格納されている。パラメータ算出部１０２は、これを用いて、衛星位置ベクトルＰ→の３成分（ｘ，ｙ，ｚ）それぞれを時刻を変数とした２次ないし２次以上のＮ次多項式で表現する（多項式の係数決定に最小二乗法を使う）。そして、それら３つのＮ次多項式に観測時刻ｔを与える。このようにすることで、観測時刻ｔにおける衛星位置ベクトルＰ→が求まる。衛星速度ベクトルＶ→については上記３つのＮ次多項式をそれぞれ時刻で微分した式（３つのＮ－１次多項式）を用意しておく。そして、それらＮ－１次多項式に観測時刻ｔを与える。このようにすることで、観測時刻ｔにおける衛星速度ベクトルＶ→が求まる。
【0043】
ステップＳ１０５－３
パラメータ算出部１０２は、視線ベクトルを算出する。以下では、視線ベクトルをＬの上に矢印（→）が付されている表記とＬの横に矢印（→）が付されている表記のいずれかにより示す。
ここでは、ステップＳ１０４で求めた観測点の３次元座標値（Ｘ，Ｙ，Ｚ）をＴの上に矢印（→）が付されている表記とＴの横に矢印（→）が付されている表記のいずれかにより示す。
パラメータ算出部１０２は、視線ベクトルＬ→、すなわち、人工衛星５００から観測点へ向かうベクトルを以下の式（３）により求める。つまり、パラメータ算出部１０２は、衛星位置ベクトルＰ→から観測点Ｔ→を減じて視線ベクトルＬ→を得る。
パラメータ算出部１０２は、視線ベクトルを算出する。以下では、視線ベクトルをＬの上に矢印（→）が付されている表記とＬの横に矢印（→）が付されている表記のいずれかにより示す。
ここでは、ステップＳ１０４で求めた観測点の３次元座標値（Ｘ，Ｙ，Ｚ）をＴの上に矢印（→）が付されている表記とＴの横に矢印（→）が付されている表記のいずれかにより示す。
パラメータ算出部１０２は、視線ベクトルＬ→、すなわち、人工衛星５００から観測点へ向かうベクトルを以下の式（３）により求める。つまり、パラメータ算出部１０２は、衛星位置ベクトルＰ→から観測点Ｔ→を減じて視線ベクトルＬ→を得る。
【0044】
【数2】
【0045】
ステップＳ１０５－４
パラメータ算出部１０２は、ドップラ角θｄを算出する。ドップラ角θｄは、視線ベクトルＬ→と衛星速度ベクトルＶ→とのなす角である。
具体的には、パラメータ算出部１０２は、以下の式（４）によりドップラ角θｄを算出する。
パラメータ算出部１０２は、ドップラ角θｄを算出する。ドップラ角θｄは、視線ベクトルＬ→と衛星速度ベクトルＶ→とのなす角である。
具体的には、パラメータ算出部１０２は、以下の式（４）によりドップラ角θｄを算出する。
【0046】
【数3】
【0047】
ステップＳ１０５－５
パラメータ算出部１０２は、ドップラ角θｄが９０°（視線ベクトルＬ→と衛星速度ベクトルＶ→との直交状態）に十分近いかどうかを判定する。
図６に示すように、９０°に対するドップラ角θｄの誤差は｜θｄ－９０°｜である。これに視線ベクトルＬ→の長さ（＝人工衛星５００から観測点までの距離）｜Ｌ→｜を乗ずれば、ゼロドップラ位置に対して人工衛星５００がどの程度ずれた位置にいるかを概算することができる。
式（５）に示すように、算出された誤差（｜θｄ－９０°｜×｜Ｌ→｜）が閾値未満であれば、パラメータ算出部１０２は、現在の観測時刻ｔ、衛星位置ベクトルＰ→及び衛星速度ベクトルＶ→をゼロドップラパラメータとして採用する。このため、処理がステップＳ１０６に進む。
パラメータ算出部１０２は、ドップラ角θｄが９０°（視線ベクトルＬ→と衛星速度ベクトルＶ→との直交状態）に十分近いかどうかを判定する。
図６に示すように、９０°に対するドップラ角θｄの誤差は｜θｄ－９０°｜である。これに視線ベクトルＬ→の長さ（＝人工衛星５００から観測点までの距離）｜Ｌ→｜を乗ずれば、ゼロドップラ位置に対して人工衛星５００がどの程度ずれた位置にいるかを概算することができる。
式（５）に示すように、算出された誤差（｜θｄ－９０°｜×｜Ｌ→｜）が閾値未満であれば、パラメータ算出部１０２は、現在の観測時刻ｔ、衛星位置ベクトルＰ→及び衛星速度ベクトルＶ→をゼロドップラパラメータとして採用する。このため、処理がステップＳ１０６に進む。
【0048】
【数4】
【0049】
一方、誤差（｜θｄ－９０°｜×｜Ｌ→｜）が閾値ε以上であれば、処理が以下のステップＳ１０５－６に進む。
閾値εは、ミリメートル単位である。閾値εは、例えば１ｍｍである。
閾値εは、ミリメートル単位である。閾値εは、例えば１ｍｍである。
【0050】
ステップＳ１０５－６
パラメータ算出部１０２は、観測時刻ｔを補正する。
つまり、パラメータ算出部１０２は、ゼロドップラ位置（真値時刻）に近づくように観測時刻ｔを補正する。
以下に観測時刻ｔの補正量の計算方法について説明する。
パラメータ算出部１０２は、観測時刻ｔを補正する。
つまり、パラメータ算出部１０２は、ゼロドップラ位置（真値時刻）に近づくように観測時刻ｔを補正する。
以下に観測時刻ｔの補正量の計算方法について説明する。
【0051】
まず、人工衛星５００が、図７に示すように、ある平面内で円運動していると仮定する。ここでの平面は、角運動量ベクトルを垂線とする平面である。ここで、角運動ベクトルはｈの上に矢印（→）が付されている表記とｈの横に矢印（→）が付されている表記のいずれかにより示す。
角運動ベクトルｈ→は、式（６）で表される。なお、式（６）の「×」はベクトルの外積を意味する。
角運動ベクトルｈ→は、式（６）で表される。なお、式（６）の「×」はベクトルの外積を意味する。
【0052】
【数5】
【0053】
観測点Ｔ→は、この平面上にはないが、観測点Ｔ→を同平面上に射影したものをＴ’→とする。なお、「Ｔ’→」は、Ｔ’の上に矢印（→）を付した表記で表することがある。
Ｔ’→は観測点Ｔ→から平面へ垂線（＝ｈ→に比例）を降ろしたものなので、式（７）で表される。なお、ｋは定数であり、これから決定する。
Ｔ’→は観測点Ｔ→から平面へ垂線（＝ｈ→に比例）を降ろしたものなので、式（７）で表される。なお、ｋは定数であり、これから決定する。
【0054】
【数6】
【0055】
式（７）とｈ→との内積をとると、式（８）が得られる。
【0056】
【数7】
【0057】
Ｔ’→は平面内のベクトルなのでｈ→と直交する。よって、式（８）の左辺はゼロとなるので、ｋが以下の式（９）のように決定される。
【0058】
【数8】
【0059】
これにより、Ｔ’→が以下の式（１０）のように定まる。
【0060】
【数9】
【0061】
平面内の各ベクトルの関係を図７に示す。
パラメータ算出部１０２は、図７において、Ｔ’→とＰ→の方向が一致するように人工衛星５００を動かす。つまり、パラメータ算出部１０２は、観測時刻ｔをＴ’→とＰ→の方向が一致する時刻に補正する。
図７において、Ｐ→とＶ→が直交する、すなわち人工衛星５００が円運動をしていると仮定すると、θとφは以下の式（１１）の関係になる。
θ＋φ＝９０° 式（１１）
パラメータ算出部１０２は、図７において、Ｔ’→とＰ→の方向が一致するように人工衛星５００を動かす。つまり、パラメータ算出部１０２は、観測時刻ｔをＴ’→とＰ→の方向が一致する時刻に補正する。
図７において、Ｐ→とＶ→が直交する、すなわち人工衛星５００が円運動をしていると仮定すると、θとφは以下の式（１１）の関係になる。
θ＋φ＝９０° 式（１１）
【0062】
しかしながら、式（１１）は厳密には成立しない。なぜならば、一般に衛星軌道は円運動では無い。しかしながら、ステップＳ１０５－５における誤差が小さくなり、演算が収束すればよく、その目的には十分に使用可能である。
【0063】
図７のθとドップラ角θｄの関係式は以下の式（１２）のようになる。
【0064】
【数10】
【0065】
式（１２）では、２つ目の等号においてＴ’→＝Ｔ→＋ｋｈ→を用いている。また、３つ目の等号において、Ｖ→とｈ→は直交（＝内積がゼロ）することを用いた（なぜならば、式（６）よりｈ→はＶ→の外積によって定義されている）。また、４つ目の等号でＴ→＝Ｐ→＋Ｌ→（式（３）より）を用いた。また、５つ目の等号でＶ→・Ｌ→＝｜Ｖ→｜｜Ｌ｜ｃｏｓθｄ（式（４）より）を用いた。
【0066】
式（１２）より、θとθｄの微小変化の関係は以下の式（１３）のように求まる（θをθｄで偏微分し、∂θ＝Δθ、∂θｄ＝Δθｄとおく）。
【0067】
【数11】
【0068】
次に、角速度ωを式（１４）で定義する。
０００００
０００００
【数12】
【0069】
図７のφの微小変化をΔφと記すと、Δφは観測時刻ｔの微小変化Δｔと角速度ωを用いて、式（１５）にように表せる。
Δφ＝ωΔｔ 式（１５）
Δφ＝ωΔｔ 式（１５）
【0070】
また、式（１３）よりΔφ＝－Δθである。よって、Δｔは、以下の式（１６）で表される。
【0071】
【数13】
式（１６）のθｄにθｄ－９０°を代入すると、θｄが９０°となるために必要な観測時刻ｔの補正量Δｔが以下の式（１７）により得られる。
【0072】
【数14】
【0073】
パラメータ算出部１０２は、観測時刻ｔをｔ＋Δｔに補正して、ステップＳ１０５－２以降の処理を行う。
【0074】
図３に戻り、ステップＳ１０６では、パラメータ算出部１０２は、スラントレンジを算出する。
スラントレンジとは、前述の通り、人工衛星５００と観測点との距離である。
パラメータ算出部１０２は、スラントレンジＲを、前述した観測点Ｔ→（楕円体高はゼロ）と衛星速度位置ベクトルＰ→を用いて以下の式（１８）で算出する。
スラントレンジとは、前述の通り、人工衛星５００と観測点との距離である。
パラメータ算出部１０２は、スラントレンジＲを、前述した観測点Ｔ→（楕円体高はゼロ）と衛星速度位置ベクトルＰ→を用いて以下の式（１８）で算出する。
【0075】
【数15】
【0076】
次に、ステップＳ１０７では、サンプル点設定部１０１が、高さ方向のサンプル点を配置する。
つまり、サンプル点設定部１０１は、ステップＳ１０２で設定したピクセル方向及びライン方向のサンプル点のそれぞれに対して、高さ方向にもＮｓａｍｐ＿ｈ個のサンプル点を配置する。Ｎｓａｍｐ＿ｈは１０～１００程度の値とする。
図８は、高さ方向にサンプル点を配置した例を示す。
つまり、サンプル点設定部１０１は、ステップＳ１０２で設定したピクセル方向及びライン方向のサンプル点のそれぞれに対して、高さ方向にもＮｓａｍｐ＿ｈ個のサンプル点を配置する。Ｎｓａｍｐ＿ｈは１０～１００程度の値とする。
図８は、高さ方向にサンプル点を配置した例を示す。
【0077】
ステップＳ１０８では、パラメータ算出部１０２が、スラントレンジ、ゼロドップラパラメータ（衛星位置、衛星速度）を用いてＥＣＲ座標値を算出する。
パラメータ算出部１０２は、地球長半径をａ（約６，３７８ｋｍ）としたとき、（ａ＋ｍｉｎ＿ｈ）から（ａ＋ｍａｘ＿ｈ）まで等間隔にＮｓａｍｐ＿ｈ回だけ地球長半径を変化させて処理を行う。
「ｍｉｎ＿ｈ」は－１０００～０［ｍｅｔｅｒ］程度である。「ｍｉｎ＿ｈ」は海面よりも低い土地の存在を想定したものである。「ｍａｘ＿ｈ」は９０００～１００００［ｍｅｔｅｒ］程度である。「ｍａｘ＿ｈ」はエベレスト山頂（８，８４９ｍ）以上の値とする。このようにすることで、地表の一番低いところから高いところまでの全てを網羅することができる。
以下の手順により観測点Ｔ～→ののＥＣＲ座標値（Ｘ，Ｙ，Ｚ）を求めることができる。観測点Ｔ～→は、前述の観測点Ｔ→とは異なり、楕円体高はゼロではない。なお、「Ｔ～→」は、Ｔの上に「～」と矢印（→）を付した表記で表することがある。
パラメータ算出部１０２は、地球長半径をａ（約６，３７８ｋｍ）としたとき、（ａ＋ｍｉｎ＿ｈ）から（ａ＋ｍａｘ＿ｈ）まで等間隔にＮｓａｍｐ＿ｈ回だけ地球長半径を変化させて処理を行う。
「ｍｉｎ＿ｈ」は－１０００～０［ｍｅｔｅｒ］程度である。「ｍｉｎ＿ｈ」は海面よりも低い土地の存在を想定したものである。「ｍａｘ＿ｈ」は９０００～１００００［ｍｅｔｅｒ］程度である。「ｍａｘ＿ｈ」はエベレスト山頂（８，８４９ｍ）以上の値とする。このようにすることで、地表の一番低いところから高いところまでの全てを網羅することができる。
以下の手順により観測点Ｔ～→ののＥＣＲ座標値（Ｘ，Ｙ，Ｚ）を求めることができる。観測点Ｔ～→は、前述の観測点Ｔ→とは異なり、楕円体高はゼロではない。なお、「Ｔ～→」は、Ｔの上に「～」と矢印（→）を付した表記で表することがある。
【0078】
以下で説明する処理は、図９に示すように、人工衛星５００の位置を中心として人工衛星５００の進行方向に直交する半径Ｒ（スラントレンジ）の円と、地球楕円体との交点を求める方法である。
【0079】
ステップＳ１０８－１
パラメータ算出部１０２は、地球長半径Ｒｅと地球短半径Ｒｐを設定する。
真の地球長半径をａ、扁平率をｆとしたとき、以降で使用する地球長半径Ｒｅと地球短半径Ｒｐを以下の式（１９）と式（２０）のように求める。ただしｉ＝０～Ｎｓａｍｐ＿ｈ－１である。また、真の地球長半径ａ、扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
パラメータ算出部１０２は、地球長半径Ｒｅと地球短半径Ｒｐを設定する。
真の地球長半径をａ、扁平率をｆとしたとき、以降で使用する地球長半径Ｒｅと地球短半径Ｒｐを以下の式（１９）と式（２０）のように求める。ただしｉ＝０～Ｎｓａｍｐ＿ｈ－１である。また、真の地球長半径ａ、扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
【0080】
【数16】
【0081】
Ｒｐ＝Ｒｅ×（１－ｆ） 式（２０）
【0082】
ステップＳ１０８－２
次に、パラメータ算出部１０２は、ドップラ角θｄを決定する。
ＡＬＯＳ－２をはじめ、現在入手可能なＳＡＲプロダクトは全てゼロドップラ結像と呼ばれる結像位置、すなわち、ＥＣＲでの衛星速度ベクトルＶ→に対して直交する方向に像を結ぶようにして画像化されている。このため、ドップラ角θｄは９０°固定とする。
次に、パラメータ算出部１０２は、ドップラ角θｄを決定する。
ＡＬＯＳ－２をはじめ、現在入手可能なＳＡＲプロダクトは全てゼロドップラ結像と呼ばれる結像位置、すなわち、ＥＣＲでの衛星速度ベクトルＶ→に対して直交する方向に像を結ぶようにして画像化されている。このため、ドップラ角θｄは９０°固定とする。
【0083】
ステップＳ１０８－３
次に、パラメータ算出部１０２は、ＥＣＲでの衛星速度ベクトルＶ→を以下の式（２１）により単位ベクトル化した単位ベクトルｕ→を用意する。単位ベクトルｕ→は、ｕの上に矢印（→）を付した表記で表することがある。
次に、パラメータ算出部１０２は、ＥＣＲでの衛星速度ベクトルＶ→を以下の式（２１）により単位ベクトル化した単位ベクトルｕ→を用意する。単位ベクトルｕ→は、ｕの上に矢印（→）を付した表記で表することがある。
【0084】
【数17】
【0085】
ステップＳ１０８－４
次に、パラメータ算出部１０２は、ＥＣＲでの衛星位置ベクトルＰ→を以下の式（２２）により単位ベクトル化した単位ベクトルｒ→を用意する。単位ベクトルｒ→は、ｒの上に矢印（→）を付した表記で表することがある。単位ベクトルｒ→は、地心方向を向く単位ベクトルである。
次に、パラメータ算出部１０２は、ＥＣＲでの衛星位置ベクトルＰ→を以下の式（２２）により単位ベクトル化した単位ベクトルｒ→を用意する。単位ベクトルｒ→は、ｒの上に矢印（→）を付した表記で表することがある。単位ベクトルｒ→は、地心方向を向く単位ベクトルである。
【0086】
【数18】
【0087】
ステップＳ１０８－５
次に、パラメータ算出部１０２は、ｕ→とｒ→を用いて、式（２３）と式（２４）の値を用意する。
次に、パラメータ算出部１０２は、ｕ→とｒ→を用いて、式（２３）と式（２４）の値を用意する。
【0088】
【数19】
【0089】
ステップＳ１０８－６
次に、パラメータ算出部１０２は、例えば、初期値としてオフナディア角θｏ＝３５°、変化量Δθｄ＝０とおく。なお、オフナディア角θｏの３５°には特に根拠は無く、２０～５０°程度の値であればどの値でもよい。
次に、パラメータ算出部１０２は、例えば、初期値としてオフナディア角θｏ＝３５°、変化量Δθｄ＝０とおく。なお、オフナディア角θｏの３５°には特に根拠は無く、２０～５０°程度の値であればどの値でもよい。
【0090】
ステップＳ１０８－７
次に、パラメータ算出部１０２は、以下の式（２５）により、ａ、ｂ、ｃ１、ｃ２の値を求める。
次に、パラメータ算出部１０２は、以下の式（２５）により、ａ、ｂ、ｃ１、ｃ２の値を求める。
【0091】
【数20】
【0092】
ステップＳ１０８－８
次に、パラメータ算出部１０２は、式（２６）により、以下のベクトルを求める。
次に、パラメータ算出部１０２は、式（２６）により、以下のベクトルを求める。
【0093】
【数21】
【0094】
ステップＳ１０８－９
次に、パラメータ算出部１０２は、ステップＳ１０８－１で設定した地球長半径Ｒｅと短半径Ｒｐを用いて、以下の式（２８）及び式（２９）のベクトルを生成する。
ここで、Ｐｘ、Ｐｙ及びＰｚは衛星位置ベクトルＶ→の各成分である。ｌｘ、ｌｙ及びｌｚは式（２７）のｌ→の各成分である。
次に、パラメータ算出部１０２は、ステップＳ１０８－１で設定した地球長半径Ｒｅと短半径Ｒｐを用いて、以下の式（２８）及び式（２９）のベクトルを生成する。
ここで、Ｐｘ、Ｐｙ及びＰｚは衛星位置ベクトルＶ→の各成分である。ｌｘ、ｌｙ及びｌｚは式（２７）のｌ→の各成分である。
【0095】
【数22】
【0096】
ステップＳ１０８－１０
次に、パラメータ算出部１０２は、式（３０）により、スラントレンジＲ～を仮算出する。式（３０）では、スラントレンジＲ～をＲの上に「～」が付されている表記で表している。
次に、パラメータ算出部１０２は、式（３０）により、スラントレンジＲ～を仮算出する。式（３０）では、スラントレンジＲ～をＲの上に「～」が付されている表記で表している。
【0097】
【数23】
【0098】
ステップＳ１０８－１１
次に、パラメータ算出部１０２は、ステップＳ１０６で求めたスラントレンジＲと、式（３０）により求めたＲ～との差の絶対値｜Ｒ－Ｒ～｜が十分に小さければ（例えば１ｍｍ）、ステップＳ１０８－１２へ移行する。
｜Ｒ－Ｒ～｜が十分に小さくない場合（例えば、１ｍｍより大きい場合）は、パラメータ算出部１０２は、以下の式（３１）により、オフナディア角θｏを修正し、ステップＳ１０８－７へ戻る。
θｏ＝θｏ＋Δθｏ 式（３１）
ここで、Δθｏは以下の式（３２）によりに求める。
次に、パラメータ算出部１０２は、ステップＳ１０６で求めたスラントレンジＲと、式（３０）により求めたＲ～との差の絶対値｜Ｒ－Ｒ～｜が十分に小さければ（例えば１ｍｍ）、ステップＳ１０８－１２へ移行する。
｜Ｒ－Ｒ～｜が十分に小さくない場合（例えば、１ｍｍより大きい場合）は、パラメータ算出部１０２は、以下の式（３１）により、オフナディア角θｏを修正し、ステップＳ１０８－７へ戻る。
θｏ＝θｏ＋Δθｏ 式（３１）
ここで、Δθｏは以下の式（３２）によりに求める。
【0099】
【数24】
【0100】
なお、ステップＳ１０８－７からステップＳ１０８－１１の繰り返しは数回で収束する。
【0101】
式（３２）の導出には、余弦定理を使用する。図１０と余弦定理より、式（３３）が得られる。
【0102】
【数25】
【0103】
式（３３）の「ａ～」（上記ではａの上に「～」）は地球半径であるが、以下の式（３４）の式展開では定数として扱う（＝すぐ下の偏微分でゼロになる）ので具体的に求める必要は無い。
式（３３）をオフナディア角θｏで偏微分する。このとき、「｜Ｐ→｜２」と「ａ～」は定数として扱い、スラントレンジＲはオフナディア角θｏに依存するものとする。
式（３３）をオフナディア角θｏで偏微分する。このとき、「｜Ｐ→｜２」と「ａ～」は定数として扱い、スラントレンジＲはオフナディア角θｏに依存するものとする。
【0104】
【数26】
【0105】
式（３４）を整理すると、式（３５）が得られる。
【0106】
【数27】
【0107】
∂θｏをΔθｏと記し、∂ＲをΔＲ＝Ｒ－Ｒ～としたものが式（３２）である。
【0108】
ステップＳ１０８－１２
最後に、パラメータ算出部１０２は、観測点「Ｔ～→」を式（３６）により求める。
最後に、パラメータ算出部１０２は、観測点「Ｔ～→」を式（３６）により求める。
【0109】
【数28】
【0110】
図３に戻り、ステップＳ１０９では、パラメータ算出部１０２は、ＥＣＲ座標値から緯度、経度及び楕円体高への変換を行う。
パラメータ算出部１０２は、ステップＳ１０８で求めた観測点「Ｔ～→」のＥＣＲ座標値（Ｘ，Ｙ，Ｚ）から緯度、経度及び楕円体高への変換は、参考文献１のＰ１７～Ｐ１８に記載されている「繰り返し計算のない方法」を使用する。
具体的には、パラメータ算出部１０２は、式（３７）に従い変換を行う。なお、以下、ＥＣＲ座標値をＸＹＺ、緯度をＢ、経度をＬ、楕円体高をＨｅ、地球長半径をａ、扁平率をｆと記す。ここでも地球長半径ａ、扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
パラメータ算出部１０２は、ステップＳ１０８で求めた観測点「Ｔ～→」のＥＣＲ座標値（Ｘ，Ｙ，Ｚ）から緯度、経度及び楕円体高への変換は、参考文献１のＰ１７～Ｐ１８に記載されている「繰り返し計算のない方法」を使用する。
具体的には、パラメータ算出部１０２は、式（３７）に従い変換を行う。なお、以下、ＥＣＲ座標値をＸＹＺ、緯度をＢ、経度をＬ、楕円体高をＨｅ、地球長半径をａ、扁平率をｆと記す。ここでも地球長半径ａ、扁平率ｆとしては、例えばＧＲＳ８０楕円体やＷＧＳ８４楕円体などの値を使用する。
【0111】
【数29】
【0112】
次に、ステップＳ１１０において、ＲＰＣ算出部１０３が、最小二乗法によりＲＰＣを算出する。
ＲＰＣ算出部１０３は、ステップＳ１０９までで得られた「ピクセル番号、ライン番号、緯度、経度、楕円体高」のデータセット（Ｎｓａｍｐ＿ｐｘ×Ｎｓａｍｐ＿ｌｎ×Ｎｓａｍｐ＿ｈ個のデータの組）を用いて、最小二乗法により有理多項式の係数（つまりＲＰＣ）を算出する。
アルゴリズムは参考文献２の式（４ａ）～（９）のｄｉｒｅｃｔｉｏｎ ｓｏｌｕｔｉｏｎ（ｎｏｎ－ｉｔｅｒａｔｉｖｅ ｓｏｌｕｔｉｏｎによる。
参考文献２：Ｔａｏ ｅｔ ａｌ． “Ａ Ｃｏｍｐｒｅｈｅｎｓｉｖｅ Ｓｔｕｄｙ ｏｆ ｔｈｅ Ｒａｔｉｏｎａｌ Ｆｕｎｃｔｉｏｎ Ｍｏｄｅｌ ｆｏｒ Ｐｈｏｔｏｇｒａｍｍｅｔｒｉｃ Ｐｒｏｃｅｓｓｉｎｇ”， ＰＨＴＯＧＲＡＭＭＥＴＲＩＣ ＥＮＧＩＮＥＥＲＩＮＧ ＆ ＲＥＭＯＴＥ ＳＥＮＳＩＮＧ， ２００１．
ＲＰＣ算出部１０３は、ステップＳ１０９までで得られた「ピクセル番号、ライン番号、緯度、経度、楕円体高」のデータセット（Ｎｓａｍｐ＿ｐｘ×Ｎｓａｍｐ＿ｌｎ×Ｎｓａｍｐ＿ｈ個のデータの組）を用いて、最小二乗法により有理多項式の係数（つまりＲＰＣ）を算出する。
アルゴリズムは参考文献２の式（４ａ）～（９）のｄｉｒｅｃｔｉｏｎ ｓｏｌｕｔｉｏｎ（ｎｏｎ－ｉｔｅｒａｔｉｖｅ ｓｏｌｕｔｉｏｎによる。
参考文献２：Ｔａｏ ｅｔ ａｌ． “Ａ Ｃｏｍｐｒｅｈｅｎｓｉｖｅ Ｓｔｕｄｙ ｏｆ ｔｈｅ Ｒａｔｉｏｎａｌ Ｆｕｎｃｔｉｏｎ Ｍｏｄｅｌ ｆｏｒ Ｐｈｏｔｏｇｒａｍｍｅｔｒｉｃ Ｐｒｏｃｅｓｓｉｎｇ”， ＰＨＴＯＧＲＡＭＭＥＴＲＩＣ ＥＮＧＩＮＥＥＲＩＮＧ ＆ ＲＥＭＯＴＥ ＳＥＮＳＩＮＧ， ２００１．
【0113】
オルソ補正部１０４は、以上により算出されたＲＰＣとＤＥＭ１６０を用いて、撮像画像１５０に対してオルソ補正を行う。
そして、オルソ補正部１０４は、オルソ補正済みの撮像画像１５０である補正済画像１７０を出力する。
そして、オルソ補正部１０４は、オルソ補正済みの撮像画像１５０である補正済画像１７０を出力する。
【0114】
＊＊＊実施の形態の効果の説明＊＊＊
本実施の形態では、ゼロドップラが成立する時刻と、当該時刻での人工衛星５００の位置及び速度とを用いてＲＰＣを算出する。そして、算出されたＲＰＣを用いてオルソ補正が行われる。このため、特許文献１の技術で必要であったマッチング処理が不要である。
また、本実施の形態では、ゼロドップラ位置と観測時刻ｔにおける人工衛星５００の位置との誤差と比較する閾値εをミリメートル単位としている。このため、ＤＥＭの精度による制約を受けない。
以上より、本実施の形態によれば、高速かつ高精度なオルソ補正が可能になる。
本実施の形態では、ゼロドップラが成立する時刻と、当該時刻での人工衛星５００の位置及び速度とを用いてＲＰＣを算出する。そして、算出されたＲＰＣを用いてオルソ補正が行われる。このため、特許文献１の技術で必要であったマッチング処理が不要である。
また、本実施の形態では、ゼロドップラ位置と観測時刻ｔにおける人工衛星５００の位置との誤差と比較する閾値εをミリメートル単位としている。このため、ＤＥＭの精度による制約を受けない。
以上より、本実施の形態によれば、高速かつ高精度なオルソ補正が可能になる。
【0115】
実施の形態２．
本実施の形態では、主に実施の形態１との差異を説明する。
なお、以下で説明していない事項は、実施の形態１と同様である。
本実施の形態では、主に実施の形態１との差異を説明する。
なお、以下で説明していない事項は、実施の形態１と同様である。
【0116】
【0117】
本実施の形態では、オルソ補正部１０４は、加工画像２５０にオルソ補正を行う。加工画像２５０は撮像画像１５０を加工して得られた画像である。加工画像２５０は、例えば、撮像画像１５０の被写体を任意の分類区分（例えば、植生の有無、災害の有無による分類区分）で分類して得られた画像である。また、加工画像２５０は、例えば、撮像画像１５０にスペックルノイズを抑制するフィルタ処理を行って得られた画像である。
実施の形態１では、オルソ補正部１０４は撮像画像１５０に対してオルソ補正を行っている。
本実施の形態では、撮像画像１５０から算出されたＲＰＣを用いて、加工画像２５０にオルソ補正を行っている。
実施の形態１では、オルソ補正部１０４は撮像画像１５０に対してオルソ補正を行っている。
本実施の形態では、撮像画像１５０から算出されたＲＰＣを用いて、加工画像２５０にオルソ補正を行っている。
【0118】
加工画像２５０も、フォアショートニングのために、そのままでは地図と重ね合わせることができない。しかし、オルソ補正部１０４が加工画像２５０にオルソ補正を行うことで、オルソ補正後の加工画像２５０である補正済画像２７０を地図と重ね合わせることができる。
【0119】
以上、実施の形態１及び２を説明したが、これら２つの実施の形態を組み合わせて実施しても構わない。
あるいは、これら２つの実施の形態のうち、１つを部分的に実施しても構わない。
あるいは、これら２つの実施の形態を部分的に組み合わせて実施しても構わない。
また、これら２つの実施の形態に記載された構成及び手順を必要に応じて変更してもよい。
あるいは、これら２つの実施の形態のうち、１つを部分的に実施しても構わない。
あるいは、これら２つの実施の形態を部分的に組み合わせて実施しても構わない。
また、これら２つの実施の形態に記載された構成及び手順を必要に応じて変更してもよい。
【0120】
＊＊＊ハードウェア構成の補足説明＊＊＊
ここで、ハードウェア構成の補足説明を行う。
プロセッサ９０１は、プロセッシングを行うＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）である。
プロセッサ９０１は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）等である。
主記憶装置９０２は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）である。
補助記憶装置９０３は、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）等である。
通信装置９０４は、データの通信処理を実行する電子回路である。
通信装置９０４は、例えば、通信チップ又はＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。
ここで、ハードウェア構成の補足説明を行う。
プロセッサ９０１は、プロセッシングを行うＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）である。
プロセッサ９０１は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）等である。
主記憶装置９０２は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）である。
補助記憶装置９０３は、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）等である。
通信装置９０４は、データの通信処理を実行する電子回路である。
通信装置９０４は、例えば、通信チップ又はＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。
【0121】
また、補助記憶装置９０３には、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）も記憶されている。
そして、ＯＳの少なくとも一部がプロセッサ９０１により実行される。
また、プロセッサ９０１はＯＳの少なくとも一部を実行しながら、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムを実行する。
ＯＳの実行により、タスク管理、メモリ管理、ファイル管理、通信制御等が行われる。
そして、ＯＳの少なくとも一部がプロセッサ９０１により実行される。
また、プロセッサ９０１はＯＳの少なくとも一部を実行しながら、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムを実行する。
ＯＳの実行により、タスク管理、メモリ管理、ファイル管理、通信制御等が行われる。
【0122】
また、サンプル点設定部１０１、パラメータ算出部１０２等の処理の結果を示す情報、データ、信号値及び変数値の少なくともいずれかが、主記憶装置９０２、補助記憶装置９０３、プロセッサ９０１内のレジスタ及びキャッシュメモリの少なくともいずれかに記憶される。
また、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムは、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ブルーレイ（登録商標）ディスク、ＤＶＤ等の可搬記録媒体に格納されていてもよい。そして、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムが格納された可搬記録媒体を流通させてもよい。
また、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムは、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ブルーレイ（登録商標）ディスク、ＤＶＤ等の可搬記録媒体に格納されていてもよい。そして、サンプル点設定部１０１、パラメータ算出部１０２等の機能を実現するプログラムが格納された可搬記録媒体を流通させてもよい。
【0123】
また、サンプル点設定部１０１、パラメータ算出部１０２等の少なくともいずれかの「部」を、「回路」又は「工程」又は「手順」又は「処理」又は「サーキットリー」に読み替えてもよい。
また、画像処理装置１００は、処理回路により実現されてもよい。処理回路は、例えば、ロジックＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＧＡ（Ｇａｔｅ Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）である。
この場合は、サンプル点設定部１０１、パラメータ算出部１０２等は、それぞれ処理回路の一部として実現される。
なお、本明細書では、プロセッサと処理回路との上位概念を、「プロセッシングサーキットリー」という。
つまり、プロセッサと処理回路とは、それぞれ「プロセッシングサーキットリー」の具体例である。
また、画像処理装置１００は、処理回路により実現されてもよい。処理回路は、例えば、ロジックＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＧＡ（Ｇａｔｅ Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）である。
この場合は、サンプル点設定部１０１、パラメータ算出部１０２等は、それぞれ処理回路の一部として実現される。
なお、本明細書では、プロセッサと処理回路との上位概念を、「プロセッシングサーキットリー」という。
つまり、プロセッサと処理回路とは、それぞれ「プロセッシングサーキットリー」の具体例である。
【0124】
最後に、本開示の諸態様を付記としてまとめて記載する。
（付記１）
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定部と、
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する画像処理装置。
（付記２）
前記パラメータ算出部は、
前記ゼロドップラが成立すると推定される時刻を成立推定時刻として設定し、前記成立推定時刻における前記飛行体の位置及び速度を算出し、前記成立推定時刻における前記進行方向ベクトルと前記視点ベクトルとのなす角であるドップラ角を算出し、前記ドップラ角に基づき、前記進行方向ベクトルと前記視点ベクトルとが直交状態にある場合の前記飛行体の位置と前記成立推定時刻における前記飛行体の位置との誤差を算出し、算出した前記誤差を閾値と比較し、
前記誤差が前記閾値未満であれば、前記成立推定時刻と、前記成立推定時刻における前記飛行体の位置及び速度とを、前記ゼロドップラパラメータとして採用し、
前記誤差が前記閾値以上であれば、前記成立推定時刻を補正し、補正後の成立推定時刻について前記飛行体の位置及び速度を算出し、前記ドップラ角を算出し、前記誤差を算出し、前記誤差を前記閾値と比較する動作を、前記誤差が前記閾値未満になるまで繰り返す付記１に記載の画像処理装置。
（付記３）
前記パラメータ算出部は、
前記閾値としてミリメートル単位の閾値を用いる付記１又は２に記載の画像処理装置。（付記４）
前記パラメータ算出部は、
サンプル点ごとに、前記飛行体とサンプル点との距離であるスラントレンジを算出し、
サンプル点ごとに、前記スラントレンジと前記ゼロドップラパラメータとを用いて、ＥＣＲ（Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｒｏｔａｔｉｎｇ）座標値を算出し、
サンプル点ごとに、ＥＣＲ座標値を緯度、経度及び楕円体高へ変換し、
前記ＲＰＣ算出部は、
サンプル点ごとに、変換により得られた緯度、経度及び楕円体高を用いて、ＲＰＣを算出する付記１～３のいずれかに記載の画像処理装置。
（付記５）
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像にオルソ補正を行うオルソ補正部を有する付記１に記載の画像処理装置。
（付記６）
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像を加工して得られた加工画像にオルソ補正を行うオルソ補正部を有する付記１～５のいずれかに記載の画像処理装置。
（付記７）
前記サンプル点設定部は、
前記撮像画像である合成開口レーダ画像に前記複数のサンプル点を設定する付記１～６のいずれかに記載の画像処理装置。
（付記８）
コンピュータが、飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
前記コンピュータが、サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
前記コンピュータが、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とを有する画像処理方法。
（付記９）
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とをコンピュータに実行させる画像処理プログラム。
（付記１）
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定部と、
サンプル点ごとに、前記飛行体の進行方向ベクトルと前記飛行体からサンプル点への視点ベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出部と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出部とを有する画像処理装置。
（付記２）
前記パラメータ算出部は、
前記ゼロドップラが成立すると推定される時刻を成立推定時刻として設定し、前記成立推定時刻における前記飛行体の位置及び速度を算出し、前記成立推定時刻における前記進行方向ベクトルと前記視点ベクトルとのなす角であるドップラ角を算出し、前記ドップラ角に基づき、前記進行方向ベクトルと前記視点ベクトルとが直交状態にある場合の前記飛行体の位置と前記成立推定時刻における前記飛行体の位置との誤差を算出し、算出した前記誤差を閾値と比較し、
前記誤差が前記閾値未満であれば、前記成立推定時刻と、前記成立推定時刻における前記飛行体の位置及び速度とを、前記ゼロドップラパラメータとして採用し、
前記誤差が前記閾値以上であれば、前記成立推定時刻を補正し、補正後の成立推定時刻について前記飛行体の位置及び速度を算出し、前記ドップラ角を算出し、前記誤差を算出し、前記誤差を前記閾値と比較する動作を、前記誤差が前記閾値未満になるまで繰り返す付記１に記載の画像処理装置。
（付記３）
前記パラメータ算出部は、
前記閾値としてミリメートル単位の閾値を用いる付記１又は２に記載の画像処理装置。（付記４）
前記パラメータ算出部は、
サンプル点ごとに、前記飛行体とサンプル点との距離であるスラントレンジを算出し、
サンプル点ごとに、前記スラントレンジと前記ゼロドップラパラメータとを用いて、ＥＣＲ（Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｒｏｔａｔｉｎｇ）座標値を算出し、
サンプル点ごとに、ＥＣＲ座標値を緯度、経度及び楕円体高へ変換し、
前記ＲＰＣ算出部は、
サンプル点ごとに、変換により得られた緯度、経度及び楕円体高を用いて、ＲＰＣを算出する付記１～３のいずれかに記載の画像処理装置。
（付記５）
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像にオルソ補正を行うオルソ補正部を有する付記１に記載の画像処理装置。
（付記６）
前記画像処理装置は、更に、
ＲＰＣを用いて、前記撮像画像を加工して得られた加工画像にオルソ補正を行うオルソ補正部を有する付記１～５のいずれかに記載の画像処理装置。
（付記７）
前記サンプル点設定部は、
前記撮像画像である合成開口レーダ画像に前記複数のサンプル点を設定する付記１～６のいずれかに記載の画像処理装置。
（付記８）
コンピュータが、飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
前記コンピュータが、サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
前記コンピュータが、サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とを有する画像処理方法。
（付記９）
飛行体により上空から撮像された、オルソ補正が行われていない撮像画像に複数の仮想点を複数のサンプル点として設定するサンプル点設定処理と、
サンプル点ごとに、前記飛行体の進行方向のベクトルと前記飛行体からサンプル点へのベクトルとが直交するゼロドップラが成立する時刻と、前記時刻での前記飛行体の位置及び速度とを、ゼロドップラパラメータとして算出するパラメータ算出処理と、
サンプル点ごとのゼロドップラパラメータを用いて、ＲＰＣ（Ｒａｔｉｏｎａｌ Ｐｏｌｙｎｏｍｉａｌ Ｃｏｅｆｆｉｃｉｅｎｔｓ）を算出するＲＰＣ算出処理とをコンピュータに実行させる画像処理プログラム。
【符号の説明】
【0125】
１００ 画像処理装置、１０１ サンプル点設定部、１０２ パラメータ算出部、１０３ ＲＰＣ算出部、１０４ オルソ補正部、１５０ 撮像画像、１６０ ＤＥＭ、１７０ 補正済画像、２５０ 加工画像、２７０ 補正済画像、５００ 人工衛星、９０１ プロセッサ、９０２ 主記憶装置、９０３ 補助記憶装置、９０４ 通信装置。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】