JP6421395B2

JP6421395B2 - Ｓａｒ図からの立体地形図形成方法

Info

Publication number: JP6421395B2
Application number: JP2014181677A
Authority: JP
Inventors: 純平上本; 清峰浦塚; 博熊谷; 孝文青木; 伊藤　康一; 康一伊藤; 酒井　修二; 修二酒井; 大樹丸木
Original assignee: Tohoku University NUC; National Institute of Information and Communications Technology
Current assignee: Tohoku University NUC; National Institute of Information and Communications Technology
Priority date: 2014-09-05
Filing date: 2014-09-05
Publication date: 2018-11-14
Anticipated expiration: 2034-09-05
Also published as: JP2016057092A

Description

この発明は、合成開口レーダ（ＳＡＲ：Synthetic Aperture Radar）によって得られる少なくとも２枚の異なる航路のＳＡＲ図から立体地形図を作成するためのＳＡＲ図からの立体地形図形成方法に関する。但し、ここで言うＳＡＲ図は、ＳＡＲ画像を含むものとする。

＜１．ＳＡＲについて＞
地球温暖化や大規模自然災害のような地球規模の問題に対処するために、リモートセンシング技術を用いた地球観測が行われている。人口衛星や航空機を用いたリモートセンシングでは、カメラのような受動型センサ、あるいは、レーダのような能動型センサが用いられる。光学センサであるカメラは、簡便に観測を行うことができるが、受動型センサであるため夜間に使用することができない大きな欠点がある。一方で、レーダは、能動型センサであるため昼夜問わずに観測することができたり、マイクロ波の透過特性から雲、霧、雨などの天候の影響を受けることなく観測できたりする利点がある。光学センサに比べて解像度が低い問題があるが、合成開口レーダ（ＳＡＲ）を用いることで、光学センサと同等以上の解像度で計測できる。本発明では、その有用性から最も重要なセンサとして利用されているＳＡＲに着目する。

ＳＡＲは、電磁波を用いて地表面の画像を生成するイメージングレーダの一種である。電磁波の位相を利用した合成開口技術により、他のイメージングレーダと比べて、ＳＡＲは高い空間分解能を持つ。人工衛星や航空機などの移動物体（プラットフォーム）にＳＡＲを搭載して、地表面などを観測する。プラットフォームは、移動しながらマイクロ波を照射し、地表面で後方散乱された反射波の強度、および、位相を測定する。

受信信号の振幅と位相情報は、受信波の特徴を表現しており、これらを用いて、複数のデータを干渉させることができる。この干渉を利用して地形変動や標高を計測する干渉ＳＡＲ(ＩｎＳＡＲ：Interferometric ＳＡＲ）が提案されている。ＩｎＳＡＲでは、ＳＡＲを用いて観測された同一地域の二つのデータの干渉から計算した位相差に基づいて、地形の標高を観測する。位相差を用いているため、生成される干渉縞の位相は、２πの不定性を持つ。また、地表が平らな場合でも軌道縞と呼ばれる干渉縞が生じる。地形図を作成するためには、軌道縞を除去する必要があるが、軌道情報の不確かさから、完全に除去することが困難であったり、除去できたとしても計測精度が低下してしまったりする場合がある。

地表の３次元計測データは、数値標高モデル（ＤＥＭ：Digital Elevation Model）が国土地理院から基盤地図情報の一つとして提供されている。これからも分かるように、この３次元計測データは、近年発達している地理情報システム（ＧＩＳ：Geographic Information System）においてその利活用が促進されているデータであり、洪水等のハザードマップの作成、都市計画、地形変化の把握等の分野で利用されている。３次元計測の有効な手段の一つとしては、飛翔体搭載センサ等によるリモートセンシングデータを用いた方法が挙げられる。これは、比較的安価であること、及び広範囲でデータを取得する能力を有することからである。その中でも、アクティブレーダであるＳＡＲは、天候、昼夜を問わず地表の観測が可能であるため、３次元計測手法の１つとして期待されている。

ＳＡＲデータを用いた３次元計測法としては、クリノメトリー（Clinometry）、ステレオスコピー（Stereoscopy、レーダ測量法（Radargrammetry）とも呼ばれる）、インターフェロメトリー（Interferometry、InSAR）、ポラリメトリー（Polarimetry）が代表的な手法として挙げられる。これらの中で最も良く用いられていると考えられる手法はＩｎＳＡＲである。しかし、ＩｎＳＡＲは位相差を用いる事により、高精度に地表面高度を計測できる長所を有するものの、適用可能なデータに厳しい観測制約条件がある事、及び手法原理上、計測結果に生得的な位相（２π）の不確定性を持つ事が短所として挙げられる。

一方、レーダ測量法は、写真測量法（Photogrammetry）をレーダに適用した手法でＩｎＳＡＲと同様に複数の図や画像を用いる手法である。計測精度がＳＡＲの画像分解能スケール程度であるため、ＳＡＲ自体の画像分解能の低さに起因して近年まではＩｎＳＡＲ手法の方が頻繁に用いられる傾向にあった。しかしながら、ＩｎＳＡＲと比較した場合、観測条件が厳しくなく、また計測結果に位相の不確定性が無い等の長所を有していることから、近年の高分解能ＳＡＲ衛星の運用開始に伴い、レーダ測量法の手法が再注目されてきている。

上記の写真測量法で本発明に近いものは、例えば特許文献１（特開２００２−２４３４４４号公報）に記載された、ステレオ計測法と呼ばれるものである。これは、通常の光学立体写真用の２枚の写真やその画像データから、被写体の立体モデルを作成することができるもので、２つの視点を用いるものである。例えば左右に離れた地点から撮像した画像における被写体のずれがカメラからの距離に依存することを用いたものである。具体的には、例えば計測対象をステレオカメラで撮像するもので、同じ対象を、視点、つまり撮像位置、を左右に移動して撮像したものである。これによって得られた画像を用いて三角測量を行う。例えば、これらの画像を重ね合わせることで、視差、つまり同じ部分の左右のズレ、を見出し、このズレの大きさから、上記対象の奥行き位置、つまり被写体の撮像カメラ位置までの近さ、を判定するものである。

また、近年の計算機能力の向上により２視点以上のディジタルカメラ画像等からマッチング技術を用いて３次元情報を復元する手法（多視点３次元復元技術）が急速に発展してきている。多視点３次元復元では、カメラの撮影位置により対象物体の画像上での写り方が変わることを利用し、画像間のマッチング結果とカメラの３次元幾何モデルから、対象物体の３次元形状を復元する。多視点３次元復元において重要となる基礎理論や要素技術としては、カメラの３次元幾何モデルに関する理論、カメラの校正やパラメータ推定技術、画像間のマッチング技術、及び複数画像の対応関係からの３次元座標計算に関する理論等が挙げられる。また、近年ではバンドル調整と呼ばれるカメラパラメータの最適化手法が確立し、事前にカメラの校正をしなくても撮影された画像のみから高精度な３次元復元を行うことができるようになってきている。

本発明に関連する従来技術としては以下が挙げられる。
Ａ．センサ/対象物間の関係式と画像マッチングを組み合わせた手法
ＳＡＲセンサと対象物間の距離と方向等から導出されるセンサ/対象物間の関係式を、画像マッチング手法から見出した複数画像上での対象物対応点に適用し、高度値を得る手法である(非特許文献１、２等)。ＳＡＲセンサの観測中の時々刻々の位置や速度等の情報を必要とする。また、高度値の精度向上の為に予め調査されたＧＣＰ（Ground Control Point）がしばしば用いられる。
Ｂ．ダイレクト立体レーダ測量法（Direct Stereo Radargrammetric）
これは、例えば非特許文献３で公開されたもので、ＳＡＲ画像の高精度な位置情報と予め生成された粗い精度のＤＥＭを必要とする。マスター画像のピクセル値に粗い精度のＤＥＭから与えられる高度に幅を持たせ、スレーブ画像上の対応点位置を逐次算出し、算出対応点近傍でのピクセル値に対する相関係数の大きさから最適な高度値を計算する手法である。
Ｃ．ディジタルカメラ画像等からのマッチング技術を用いた３次元復元の従来技術
これは、カメラの撮影位置と対象物体形状に依存するカメラ画像間の投影位置の違いから、対象物体の３次元形状を計算するコンピュータビジョンに関するもので、非特許文献４、５等の例がある。この手法を適用するためには、カメラパラメータと呼ばれるカメラの焦点距離や撮像素子の大きさ、撮影位置・姿勢等を含む情報を必要とする。これらの情報は、既知パターンを用いたカメラ校正や、複数画像間の対応関係からのパラメータ推定により取得する。カメラパラメータ推定では、高精度化のためにバンドル調整がしばしば用いられる。また、複数画像間の対応関係を求めるために、カメラ画像間のマッチングを行う必要がある。

カメラ画像間のマッチング技術としては、ＳＩＦＴ（Scale Invariant Feature Transform）に代表される特徴ベースのマッチングと、ＳＡＤ（Sum of Absolute Difference）やＰＯＣ（Phase-Only Correlation)に基づく手法に代表される領域ベースのマッチングがある。一般に、カメラパラメータの推定には特徴ベースのマッチングが、３次元形状の復元には領域ベースのマッチングが用いられる。

このように、コンピュータビジョンの分野では、複数枚のカメラ画像を用いて３次元計測を行う手法が数多く提案されている。ＳＡＲ画像を、従来のディジタルカメラ画像のように扱うことができれば、コンピュータビジョンの原理に基づいて、複数の航路で取得したＳＡＲ画像から３次元計測を行うことが可能である。しかし、ＳＡＲ画像とカメラ画像とでは画像生成プロセスが異なるため、カメラ画像を用いた３次元計測手法をＳＡＲ画像に適用することはできない。
＜２．ＳＡＲ画像の生成と特徴＞
ＳＡＲを用いて地表面を計測した信号から画像を生成するための原理とＳＡＲ画像の特徴について述べる。なお、本発明では、図１に示すように、ＳＡＲ観測でよく用いられる基本的な用語を用いる。レーダを搭載する人口衛星や航空機をプラットフォーム、プラットフォームの進行方向をアジマス方向、進行方向と直角でマイクロ波を照射する方向をレンジ方向と呼ぶ。さらに、レンジ方向は、マイクロ波が照射される方向に対するスラントレンジ方向と、地表面を基準としたグランドレンジ方向に区別される。また、マイクロ波が照射される領域で、アンテナに近い側をニアレンジ、遠い側をファーレンジと呼ぶ。
＜２．１ＳＡＲ画像の生成＞
次に、ＳＡＲの強度信号から画像を生成する原理について図１に沿って説明する。プラットフォームからレンジ方向にマイクロ波パルスを送信し、地表面で後方散乱された反射波を受信する。この一連の動作を、アジマス方向に移動しながら繰り返すことによって２次元平面を走査する。図２に受信信号の例を示す。この様な１次元時間関数として表される受信信号を、それぞれの送信信号からの遅れに関する時間関数として並べ替えることにより２次元画像を生成する。図３に、生成した２次元画像の例を示す。各画素の輝度値が受信信号の振幅（強度）であるため、生成された画像は、レーダ強度画像と呼ばれる。実際は、分解能を高めるために、アジマス方向、レンジ方向のそれぞれについてＳＡＲ圧縮処理を行う必要があるが、この圧縮処理の詳細については、例えば非特許文献６に記載がある。
＜２．２ＳＡＲ画像の特徴＞
次に、ＳＡＲ画像には、画像の生成プロセスに起因するいくつかの特徴として、(i)アジマス方向とレンジ方向で分解能が大きく異なる場合があること、(ii)レーダ観測特有の画像変調が生じていることが挙げられる。

一般に、ＳＡＲは、レンジ方向に比べアジマス方向の分解能が高い場合が多い。また、グランドレンジ方向では、ニアレンジ側からファーレンジ側にかけて分解能が向上していく。そのため、生成されたＳＡＲ画像は、地表面を斜め上空から見たような画像となる。実際の地表面と対応させるためには、画像全体の分解能を揃え、地表面を真上から見たような俯瞰画像に変換する必要がある。なお、この俯瞰画像への変換を地上投影変換、地上投影変換前の画像をスラントレンジ画像、地上投影変換後の画像をグランドレンジ画像と呼ぶ。

また、ＳＡＲ画像特有の画像変調は、ＳＡＲの幾何学的効果としてジオメトリック画像変調と呼ばれる。ジオメトリック画像変調には、フォアショートニング、レイオーバー、陰影効果がある。
＜フォアショートニング、レイオーバー＞
ＳＡＲは、マイクロ波を斜め下方向に照射し、地表面に散乱されて戻ってきた順に反射波を記録する。そのため、プラットフォームまでの距離が短い地点ほど、プラットフォーム側の画像上に投影される。同じ水平位置でも標高が高いほどプラットフォームまでの距離が短くなるため、高さがある物体は、本来の位置よりプラットフォーム側に投影される。この現象をフォアショートニングと呼ぶ。フォアショートニングがさらに大きくなると、図４に示すように、散乱体ＡとＢとが本来の位置関係と逆転して画像上に投影されてしまう。この現象をレイオーバーと呼ぶ。フォアショートニング、レイオーバーは、対象物の高さや、レーダの入射角の大きさに依存する。
＜陰影効果＞
図４のように、充分な高さがある物体（例えば山）によって、アンテナからのマイクロ波が遮られると、その物体の点Ｂより後ろ側の領域には、マイクロ波が照射されない。そのため、この領域からの受信信号は０となり、ＳＡＲ画像上に影のように映ってしまう。この現象は陰影効果と呼ばれ、地表法線に対するレーダ波の入射角が大きくなるにつれて、この陰影効果の影響を受ける領域が増加する。しかしこの場合、フォアショートニングが小さくなる。そのため、ＳＡＲ画像の利用目的に合わせ、適切なレーダ波の入射角を選択して観測を行う必要がある。

特開２００２−２４３４４４号公報

Toutin, T., "Radarsat-2 DSM Generation With New Hybrid, Deterministic, and Empirical Geometric Modeling Without GCP", IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, Vol. 50., 2049-2055, 2012. Capaldo, P., Crespi, M., "Fratarcangeli, F., Nascetti, A., and Pieralice, F., High-Resolution SAR Radargrammetry: "A First Application With COSMO-SkyMed SpotLight Imagery", IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, Vol 8., 1100-1104, 2011. Balz, T., Zhang, L., and Liao, M.S., "Direct stereo radargrammetric processing using massively parallel processing", ISPRS JOURNAL OF PHOTOGRAMMETRY AND REMOTE SENSING, Vol. 79, 137-146, 2013. R. Szeliski, "Computer Vision: Algorithms and Applications", Springer-Verlag New York Inc., 2010. Y. Furukawa, "Multi-view 3D Reconstruction Techniques in Computer Vision", 情報処理学会研究報告コンピュータビジョンとイメージメディア, Vol. 176, No. 12, pp. 1-9, 2011. 大内和夫：リモートセンシングのための合成開口レーダの基礎，東京電機大学出版局 (2004). Takita, K., Aoki, T., Sasaki, Y., Higuchi, T. and Kobayashi, K.: "High-accuracy subpixel image registration based on phase-only correlation", IEICE Trans. Fundamentals, Vol. E86-A, No. 8, pp. 1925-1934, 2003. 青木孝文,伊藤康一,柴原琢磨,長嶋聖, "位相限定相関法に基づく高精度マシンビジョン",IEICE Fundamentals, Rev, 2007. Takita, K., Muquit, M. A., Aoki, T. and Higuchi, T., "A Sub-Pixel Correspondence Search Technique for Computer Vision Applications", IEICE Trans. Fundamentals, Vol. E87-A, No. 8, pp. 1913-1923, 2004. Dellinger, F., Delon, J., Gousseau, Y., Michel, J. and Tupin, F.: SAR-SIFT: "A SIFT-like algorithm for applications on SAR images", Proc. Int’l Geoscience and Remote Sensing Symp., pp. 3478-3481 (2012).

従来技術の大きな課題の一つとして、画像処理に必要となるパラメータを予め取得することが困難であることが挙げられる。例えば、非特許文献３の手法では粗い精度で良いとは言え、予め３次元計測対象領域のＤＥＭが必要となる。このため、この従来技術は、大きく地形が変化した場所（例えば西之島）では適用不能と考えられる。また、非特許文献２においてはＳＡＲ観測中のプラットフォームの時々刻々の速度と位置情報が必須とされ、ＧＣＰ無しでも計算はできるものの、その場合は急激に精度が落ちる。更に非特許文献１では、ＧＣＰは必要としないものの、画像上のピクセル位置と実際の緯度／経度／高度を結び付けるための近似式（ＲＦＭ）、あるいはその近似式パラメータ（ＲＰＣｓ）が計算に必要とされるが、非特許文献２によると代表的な商用衛星ＳＡＲであるTerraSAR-X、COSMO-Skymedについては、ＲＦＭ、ＲＰＣｓが提供されていない。このため、３次元計測においては、予め必要となるパラメータの数が少なく、かつ取得が容易である手法が望まれる。

また、ディジタルカメラ画像等からのマッチング技術を用いた３次元復元手法はカメラ画像を対象としており、ＳＡＲ画像に直接利用することができない点が課題として挙げられる。多視点３次元復元をＳＡＲに応用するためには、ＳＡＲ画像における投影プロセスのモデル化、及び各基礎理論や要素技術のＳＡＲ画像に特化した形での構築や開発が必要である。

そこで、本発明では、少なくとも２枚の異なる航路のＳＡＲ図から立体地形図を作成するための、ＳＡＲ図からの立体地形図形成方法、を提案する。

本発明では、ＳＡＲ画像をカメラ画像と同様に扱うためにＳＡＲの幾何学的関係を定式化し、コンピュータビジョンの原理に基づいて３次元計測を行う手法を提案する。具体的には、ＳＡＲの画像生成プロセスをレーダ画像投影モデルとして導出し、ＳＡＲにおける内部または外部パラメータの定義を行う。これにより、三角測量に基づく３次元計測、バンドル調整によるパラメータの最適化を可能とする。提案手法は、ＩｎＳＡＲによる標高計測とは異なり、位相情報を利用しないので、位相を利用することによる問題点（２πの不定性、軌道縞の問題、プラットフォームの軌道制約）を回避しつつ、３次元計測を行うことが可能である。

このように、本発明の場合は、視差を用いるものではなく、奥行きを判定するものでもない。本発明は、ＳＡＲ図における、フォアショートニング効果等を用いて、計測対象の高さを測定するものであり、上記のステレオ計測法とは、用いる原理や、測定対象とする長さの方向が異なっている。

より詳しく説明すると、本発明のＳＡＲ図からの立体地形図形成方法は、同一地域のデータを含む少なくとも２枚の互いに異なる航路で取得されたＳＡＲ図を用いて、上記同一地域の所定地点の標高を含む座標を決定する方法であって、上記所定地点において、地表の対象点の位置のスラントレンジ面への射影における歪みを測定して上記標高を含む座標を決定するもので、以下の特徴を備えるものである。つまり、
（１）少なくとも２枚の上記ＳＡＲ図を、同一方角に向き同一縮尺のグランドレンジ面上のＳＡＲ図に変換するステップと、
（２）グランドレンジ面上の上記ＳＡＲ図から共通の領域を選択するステップと、
（３）共通の上記ＳＡＲ図の一方における上記所定地点付近の局所図を動かして、上記ＳＡＲ図の他方との最大一致を示す点を見いだすことで、上記所定地点の異なる上記ＳＡＲ図上での対応点を見いだすステップと、
（４）見いだした各対応点の座標をスラントレンジ面上の座標に変換して、スラントレンジ面上の該対応点の座標の位置から上記所定地点の標高値を見いだすステップと、
を含む。

上記所定地点のグランドレンジ面上のＳＡＲ図における位置としては、上記グランドレンジ面上の上記局所図の重心位置とすることができる。

上記所定地点は、上記同一地域を構成する複数の各点とすることができる。つまり、図面の全ての地点の標高を求めることができる。

特に、本発明は、同一地域のデータを含む少なくとも２枚の互いに異なる航路で取得されたＳＡＲ図を用いて、上記同一地域の所定地点の標高を含む座標を推定することを含むＳＡＲ図からの立体地形図形成方法であって、上記の特徴の他に、次の特徴を備えるものであってもよい。つまり、
（１）上記ＳＡＲ図の各々にグランドレンジ変換を行って各々のグランドレンジ図を得るステップと、
（２）上記各々のグランドレンジ図について、それらを重ね合わせるために必要な図間の回転角と並行移動量を推定するステップと、
（３）所定の地点に対応する上記各々のグランドレンジ図上の地点の座標を要素とする対応点群を求めるステップと、
（４）上記対応点の各々の要素について、上記グランドレンジ変換の逆変換を行ってスラントレンジ図上の対応点群に座標変換するステップと、
（５）スラントレンジ図上の上記対応点群について、プラットフォーム高度、スラントレンジ方向のビーム入射角、スラントレンジ方向、およびアジマス方向の分解能の逆数を含む上記ＳＡＲ図の内部パラメータと、回転角と並行移動量とを含む外部パラメータと、を用いて上記対応点群の各々の要素の高度を推定するステップと、
（６）上記内部パラメータまたは／および外部パラメータについて、バンドル調整に基づく再投影誤差の最小化を実施し、該再投影誤差が最小となる最適内部パラメータまたは／および最適外部パラメータを求めるステップと、
（７）上記最適内部パラメータまたは／および最適外部パラメータを用いて、高度推定値を得るステップと、
を含む。

上記の（２）上記各々のグランドレンジ図について、それらを重ね合わせるために必要な図間の回転角と並行移動量を推定するステップと（３）の所定の地点に対応する上記各々のグランドレンジ図上の地点の座標を要素とする対応点群を求めるステップでは、位相限定相関法（ＰＯＣ：Phase-Only Correlation）を適用することができる。

上記ＳＡＲ図は、特に、単偏波スラントレンジ振幅強度図である場合に、発明の利点が明確である。

上記バンドル調整における評価関数としては、種々の形式から選択することができるが、簡明な例としては、以下のものを挙げることができる。つまり、評価関数としては、少なくとも２枚の上記ＳＡＲ図のスラントレンジ面上の対応点群の各要素の座標と上記ＳＡＲの内部パラメータおよび外部パラメータを用いて対応点群の各要素の３次元座標を算出し、再度２次元スラントレンジ面に投影した際の対応点群の各要素の座標との差の２乗和とする。

本発明の場合、干渉位相を用いずにＳＡＲ強度画像から３次元高度情報を得る事ができる。その際、ＩｎＳＡＲのように厳しい観測制限は課されない。また、予めＧＣＰやＤＥＭ等の地上情報を必要としない。

ＳＡＲ観測でよく用いられる基本的な用語を示す図である。プラットフォームからレンジ方向にマイクロ波パルスを送信し、地表面で後方散乱された反射波を受信する際の、受信信号の例を示す図である。受信信号から生成した２次元画像の例を示す図である。陰影効果の例を示す図である。レーダ画像投影モデルを定義するための、３次元空間と２次元空間のそれぞれの位置を表す座標系を示す図である。スラントレンジ画像間の対応点を取得するまでのフローを示す図である。（ａ）マスタースラントレンジ画像（左）とスレーブスラントレンジ画像（右）、（ｂ）グランドレンジ変換画像、（ｃ）グランドレンジ画像上の対応点の表示例、（ｄ）グランドレンジ画像間の対応点の取得、（ｅ）スラントレンジ画像上の対応点の表示例、（ｆ）高度計測結果（最適化前）、（ｇ）高度計測結果（最適化後）実験で使用したＳＡＲ画像で、（ａ）が２０,０００×１０,７１３ピクセル、（ｂ）が２０,０００×１０,３３８ピクセルである。地上投影変換後のグランドレンジ画像と標高マップを示し、（ａ）および（ｂ）は、それぞれ航路１および航路２の地上投影変換後のグランドレンジ画像であり、（ｃ）は、画像マッチングの際に基準点を配置した領域をグランドレンジ画像から抽出したものであり、（ｄ）は、バンドル調整前のパラメータを用いて計算した標高マップであり、（ｅ）は、バンドル調整で最適化されたパラメータを用いて計算した標高マップであり、（ｆ）は、国土地理院の基盤地図情報（数値標高モデル）データ５ｍメッシュ（標高）を使用した標高マップである。実施例１を単純化した実施例２の配置を示す図である。

以下に、この発明の実施の形態を図面に基づいて詳細に説明する。

汎用性のある３次元計測手順の概略を以下に示す。

本実施例の手法は、例えば、ＳＡＲ画像として最もプリミティブなプロダクトである単偏波スラントレンジ振幅強度画像に適用可能である。本手法を適用するために必要となる画像枚数は２枚以上である。ここでは説明を簡略化するため、２枚の画像を用いる際の手法について記載する事とし、一枚の画像をマスター画像、残りの一枚をスレーブ画像と称する。
＜手順１＞グランドレンジ変換（画像マッチングの下準備）
手順２における一連の画像処理の下準備として、マスター画像、スレーブ画像の両スラントレンジ画像（図６（ａ））について地上投影変換を実施する。この変換は、グランドレンジ変換と呼ばれる一般的な変換である。このグランドレンジ変換により得られる画像は、グランドレンジ画像と呼ばれており、画像全体での分解能が揃っているため、異なるＳＡＲ画像同士での画像マッチングが容易となる（図６（ｂ））。なお、このグランドレンジ変換の際に必要となる各ピクセルの高度Ｚ_Aについてはこの段階では任意で良い。また、グランドレンジ変換画像が既に手元にある場合は、この手順はスキップしても良い。
＜手順2＞外部パラメータの推定と画像マッチング
手順１にて変換や生成されたグランドレンジ画像ペアについて、非特許文献８により公知である回転不変ＰＯＣを適用し、画像間回転角φを推定する。その後、スレーブ画像のグランドレンジ画像について回転角φだけ回転させ、非特許文献７により公知であるＰＯＣを用いて並行移動量（ｔ_x、ｔ_y）を推定する。推定した回転角φと並行移動量（ｔ_x、ｔ_y）を用いて変換したスレーブ画像とマスター画像について非特許文献９により公知であるＰＯＣに基づく対応点検索を実施し、グランドレンジ画像間の対応点を求める（図６（ｃ））。
＜手順３＞対応点のスラントレンジ変換
手順２で求められたグランドレンジ画像上の対応点のピクセル座標値について手順１と逆の手順でスラントレンジ画像上のピクセル座標値に変換する。これにより手順２で求めた対応点がスラントレンジ画像上に配される（図６（ｄ））。
＜手順４＞内部および外部パラメータを利用した対応点の初期高度計測
手順３までで得られたスラントレンジ画像上の対応点について、ＳＡＲ画像とともに提供される内部パラメータ（プラットフォーム高度Ｚ₀、スラントレンジ方向のビーム入射角θ、スラントレンジ方向、アジマス方向の分解能の逆数α_u、α_v）、及び手順２で推定した外部パラメータ（回転角φと並行移動量（ｔ_x、ｔ_y））を用いて、数７を計算する事で対応点における高度推定を実施する（図６（ｅ））。
＜手順５＞内部または／および外部パラメータの最適化
手順２で推定した外部パラメータは手順１で画像全体に渡り任意高度Ｚ_Aを持つとしてグランドレンジ変換を実施しているため、真値では無い。また、ＳＡＲ画像とともに提供されるパラメータについても誤差を含む可能性がある。これらの内部または／および外部パラメータについて非特許文献４の記載に沿ってバンドル調整に基づく再投影誤差の最小化を実施し、その再投影誤差が最小となる内部または／および外部パラメータを求める。
＜手順６＞最適化された内部および外部パラメータに基づく高度計測
手順５で求められた内部および外部パラメータを用いて、最終的な高度推定値を得る（図６（ｆ））。
＜３．ＳＡＲの３次元幾何の定式化＞
ここでは、コンピュータビジョンの原理に基づいて３次元計測を行うために、コンピュータビジョンの分野における解法に合わせてＳＡＲの３次元幾何を定式化する。

まず、３次元空間上の物体が２次元レーダ画像上にどのように投影されるかを記述する投影モデルとＳＡＲの内部パラメータを定義する。次に、異なる航路で同一の領域をＳＡＲで観測した場合のＳＡＲの外部パラメータを定義する。最後に、２枚のレーダ強度画像の対応関係とＳＡＲの３次元幾何に基づいて、地表面の３次元形状を求めるための理論式を導出する。
＜３．１レーダ画像投影モデルとＳＡＲの内部パラメータ＞
レーダ画像投影モデルを定義するために３次元空間と２次元空間のそれぞれの位置を表す座標系（図５）を導入する。２次元空間の座標系には、原点を画像左上とし、アジマス方向を水平軸ｕ、レンジ方向を垂直軸ｖとするディジタル画像座標系を用いる。３次元空間の座標系には、コンピュータビジョンの分野で用いられるカメラ座標系にならって、レーダ座標系と呼ぶ新たな座標系を導入する。レーダ座標系では、ディジタル画像座標系の原点に対応する３次元空間位置の標高０ｍの点を原点とし、アジマス方向をＹ_R軸、レンジ方向をＸ_R軸、高さをＺ_R軸とする。このとき、２節で述べたＳＡＲによる画像生成プロセスを考慮すると、レーダ座標系で表される３次元空間上の点（Ｘ_R、Ｙ_R、Ｚ_R）とその投影点（ｕ、ｖ）との間に次式で表される投影モデルが定義できる。

ここで、Ｘ₀はプラットフォームから原点までの水平距離、Ｄ_SLはプラットフォームから原点までのスラントレンジ距離、Ｚ₀はプラットフォームの高度を表す。また、α_uとα_vは、それぞれアジマス方向とスラントレンジ方向に関する分解能の逆数を表す定数である。Ｘ₀とＤ_SLは、Ｚ₀とレーダ入射角θを用いて以下の式で表すことができる。

以上より、レーダ画像投影モデルの計算に必要なパラメータは、プラットフォーム高度Ｚ₀、スラントレンジ方向のビーム入射角θ、スラントレンジ方向、アジマス方向それぞれの分解能の逆数α_u、α_vである。これらを内部パラメータと呼ぶ。これらの内部パラメータは通常ＳＡＲ画像とともに提供されるパラメータである。

数１はアジマス方向の画像投影位置とレンジ方向の画像投影位置とを表す。上記のＳＡＲ画像の生成についての説明で述べたように、ＳＡＲによるアジマス方向の画像生成は、１次元時間関数として表される受信信号を、送信信号から遅れの時間の関数として並べ替えることによって行われる。つまり、数１の第１式が示すように、ディジタル画像座標ｕは、単純にレーダ座標Ｙ_Rの定数倍で表現することができる。一方で、レンジ方向の画像投影位置は、プラットフォームから対象物までのスラントレンジ距離によって決定される。数１の第２式の第一項は、プラットフォームから（Ｘ_R、Ｙ_R、Ｚ_R）座標上の物体までの、次の数３の距離を、表している。

数１の下段の式で、第１項から第２項Ｄ_SLを引くことにより、（Ｘ_R、Ｙ_R、Ｚ_R）＝０のときに、（ｕ、ｖ）＝０となり、レーダ座標系とディジタル画像座標系の原点が一致する。
＜３．２座標系間の変換式とＳＡＲの外部パラメータ＞
異なる２つの航路１と航路２で同一の領域を観測した場合を考える。この時、それぞれの航路で定義されるレーダ座標系における３次元座標（Ｘ_R1，Ｙ_R1，Ｚ_R1）と（Ｘ_R2，Ｙ_R2，Ｚ_R2）との間で次式が成り立つ。

ここで、Ｒは回転角φの３×３回転行列、ｔは並進ベクトルである。ここで、数４が示す座標系間の変換式を記述するために必要なＲ、ｔをＳＡＲの外部パラメータと定義する。ＳＡＲの観測が行われている間、プラットフォームの高度は一定であると仮定すると、ＲはＺ_R軸回りの回転、ｔはＸ_R、Ｙ_R方向の並進のみを考慮するだけでも良い。この時、Ｒ、ｔは，以下の式で表すことができる。

＜３．３３次元計測の理論式＞
導出した投影モデル（数１）が、ある世界座標系の３次元座標（Ｘ，Ｙ，Ｚ）の、異なる２つの画像上のディジタル画像座標（ｕ，ｖ）と（ｕ´，ｖ´）とでそれぞれ成り立つことから、異なるレーダ座標系間の変換式（数４）を用いて、これらのディジタル画像座標（ｕ，ｖ）と（ｕ´，ｖ´）との間に，以下の式が成り立つ。

ここで，（ｕ，ｖ），（ｕ´，ｖ´）を対応点ペアとし，（Ｘ_R2，Ｙ_R2，Ｚ_R2）を世界座標系（Ｘ，Ｙ，Ｚ）とした。

数６の４つの方程式をＸ、Ｙ、Ｚについて解くことで、３次元座標を算出することができる。解法は、いくつか知られているが、最も簡明な例として、以下の方法を挙げる。つまり、上記４つの方程式から３つの方程式を選択することで、方程式を解析的に解く。この場合、３次元座標Ｘ、Ｙ、Ｚは以下の式で表すことができる。

以上のように導出されたＸ、Ｙ、Ｚは、（Ｘ_R2，Ｙ_R2，Ｚ_R2）を世界座標系と仮定した場合の３次元計測の理論式である。ただし、この仮定は本発明において必須ではない。つまり、世界座標とレーダ座標が異なっていたとしても、計測されたレーダ座標系における３次元座標を、ある剛体変換で世界座標系の３次元座標に変換することが可能である。
＜４．レーダ強度画像からの地表面の３次元計測＞
レーダ強度画像間の対応関係とＳＡＲの３次元幾何に基づいて３次元計測を行う手法を詳しく説明する。その手法は、（１）地上投影変換、（２）外部パラメータの推定と画像マッチング、（３）内部または／および外部パラメータの最適化、（４）３次元計測の４つの処理で構成される。図６に、スラントレンジ画像間の対応点を取得するまでのフローを示す。以下では、それぞれの処理について具体的に説明する。
＜４．１地上投影変換＞
ＳＡＲ画像は、上記の２．（ＳＡＲ画像の生成と特徴）〜２．２（ＳＡＲ画像の特徴）の欄で述べたように、画像全体で分解能が変化しており、実際の地表面の見え方と異なる。航路によって被写体の映り方の違いが大きくなるため、そのままでは画像マッチングを行うことが困難である。そこで、地上投影変換により画像全体の分解能を揃えて俯瞰画像に変換することで、画像間の大きな変形を解消する。（ｕ、ｖ）を地上投影変換前のスラントレンジ画像上の点とし、（Ｘ_R、Ｙ_R）を地上投影変換後のグランドレンジ画像上の点とする。この時、地上投影変換式は、２節で定義したレーダ画像投影モデルである数１で表すことができる。ただし、Ｚ_Rは任意の定数であるとする。本実施例ではＺ_R＝０とした。

この地上投影変換式により、航路１と航路２のスラントレンジ画像をそれぞれのグランドレンジ画像へ変換する。
＜４．２外部パラメータの推定と画像マッチング＞
グランドレンジ画像から外部パラメータの推定と画像マッチングを行う。まず、画像間の回転および平行移動を補正する。外部パラメータである回転行列Ｒと並進ベクトルｔは、Ｚ方向の回転と並進が０である。これを利用することで、２次元のグランドレンジ画像間の回転量と平行移動量から外部パラメータを求めることができる。本発明では、画像間の回転量推定に回転不変位相限定相関法（ＲＩＰＯＣ:Rotation-invariant Phase-Only Correlation）を、並行移動量の推定にＰＯＣを用いる（非特許文献７）。ＰＯＣは、並行移動のみの画像変形を仮定している手法であるため、はじめにＲＩＰＯＣを用いて画像間の回転角度を求める。また、航路２の画像を回転補正してからＰＯＣを用いて画像間の並行移動量を求める。次に、航路１のグランドレンジ画像と、回転補正した航路２のグランドレンジ画像との間で画像マッチングを行う。航路１のグランドレンジ画像上に基準点を配置し、回転補正した航路２のグランドレンジ画像上の対応点をＰＯＣに基づく対応点探索（非特許文献９）を用いて求める。得られた対応点を回転補正前の座標系に戻すことで、航路１と航路２のグランドレンジ画像上での対応点を得ることができる。最後に、グランドレンジ画像間の対応点ペアに、地上投影変換と逆の変換を行うことで、各航路のスラントレンジ画像上の対応点を求める。
＜４．３パラメータの最適化と３次元計測＞
グランドレンジ画像間の回転角度と平行移動量から求めた外部パラメータは、厳密に正しい値ではない。地上投影変換でＺ_R＝０を仮定しているので、画像間の平行移動量に、航路の違いによる平行移動量とジオメトリック画像変調による平行移動量が含まれてしまい、両者の区別がつかないことが原因である。また、内部パラメータには、計測による誤差が含まれていると考えられる。そこで、バンドル調整に基づく再投影誤差の最小化により、内部または／および外部パラメータの最適化を行う。

最適化を行う内部または／および外部パラメータをＰ、各航路におけるスラントレンジ画像上のＮ個の対応点ペアをｍ_1、i＝（ｕ₁、ｖ₁）_i、ｍ_2、i＝（ｕ₂、ｖ₂）_i、（１≦ｉ≦Ｎ）とする。このとき、次式で定義されるコスト関数Ｅ（Ｐ）の最小化を行う。

ここで、ｍ_rep、i（Ｐ）は、対応点ペアｍ_1、i、ｍ_2、iとパラメータＰを用いて、３次元計測の理論式に基づいて復元した３次元点をレーダ画像投影モデルにより再度２次元平面に投影した点の座標（ｕ_rep、i、ｖ_rep、i）を表す。また、ｍ_rep1、i（Ｐ）は、航路１のスラントレンジ画像上に再投影した点の座標を表し、ｍ_rep2、i（Ｐ）は、航路２のスラントレンジ画像上に再投影した点の座標を表す。

航路ｉを世界座標系として定義した場合の３次元計測の理論式を利用して求めたコスト関数をＥ_i（Ｐ）とする。このとき、最小化を行うコスト関数Ｅ（Ｐ）を以下のように定義する。

数９で定義されるコスト関数は非線形関数である。このコスト関数の最小化を行う方法は、多変数関数の最適化問題として既によく知られており、一例として、ここでは、非線形最小二乗アルゴリズムの１つであるレーベンベルグ−マルカート（Levenberg-Marquardt）法によりコスト関数の最小化を行う。そして、コスト関数Ｅ（Ｐ）が最小となるときのパラメータＰを最適化されたパラメータとして用いる。以上により求めた内部および外部パラメータとスラントレンジ画像間の対応点ペアを用いて、数７に基づいて３次元計測を行う。

上記実施例の特殊な例として、図９に示すように、２つの航路１と航路２が直交している場合の簡単なケースについて、３次元計測を行う例を示す。これは、数５の余弦項がゼロとなり、正弦項が−１となる場合である。また、簡単のため、ここでは航路１、航路２のそれぞれで観測した場合のディジタル画像（ｕ´，ｖ´）、（ｕ，ｖ）における対象物の位置はそれぞれ同定されているものとする。

この場合の、ディジタル画像と３次元空間の変換式は以下の通りである。

航路１、航路２に対する上記変換式は、それぞれ以下の数１１、数１２の通りである。

図９からＲ１系とＲ２系の関係は、この例ではφ＝−９０°に相当し、（Ｘ_R2，Ｙ_R2，Ｚ_R2）≡（Ｘ，Ｙ，Ｚ）とすると、次式となる。

数１１、１２、１３から次の式を得る。

数１４から対象物の次の３次元座標を得ることができる。

航路１が航路２と直交する様にすることで、ディジタル画像（ｕ，ｖ）、（ｕ´，ｖ´）における対象物の位置を同定することが容易になり、３次元座標を容易に求めることができる、という利点がある。

しかし、本実施例では、航路１、２が直交から２°ずれた場合、標高が１０％程度ずれる。この様に、この角度のずれが標高誤差に与える影響が比較的大きいため、上記のバンドル調整は重要である。

＜５．精度評価実験＞
本発明で提案する２枚のＳＡＲ画像から３次元計測を行う手法の計測精度を国土地理院より公開されている数値標高モデルと比較することで評価する。
＜５．１データセット＞
本実験で用いるＳＡＲ画像は、情報通信研究機構が開発した航空機ＳＡＲであるＰｉ−ＳＡＲ２により取得されたものである。観測場所は茨城県東茨城郡城里町であり、観測領域はおよそ５ｋｍ×５ｋｍである。本実験で使用したＳＡＲ画像を図７に示す。ＳＡＲ画像の大きさは、図７（ａ）が２０,０００×１０,７１３ピクセル、図７（ｂ）が２０,０００×１０,３３８ピクセルである。これらのＳＡＲ画像を取得したときのＰｉ−ＳＡＲ２の撮影諸元を表１に示す。

＜５．２実験方法＞
まず、本発明を用いて２枚のＳＡＲ画像から地表面の３次元計測を行う。そして、復元した３次元点群と、真値であるＤＥＭの３次元メッシュモデルをＩＣＰ（Iterative Closest Point）により位置合わせをし、標高の残差で評価する。国土地理院によって公開されているＤＥＭは、航空レーザー測量によって計測したデータから、建物および橋などの人工構造物や樹木などの植生をフィルタリング処理によって除去し、５ｍ間隔に標高を内挿補間して求めた数値標高モデルデータである。本発明によって求められる３次元点群が１ｍ間隔であるのに対し、ＤＥＭは５ｍ間隔の点群である。そのため、３次元点群であるＤＥＭをドロネー三角形分割により３次元メッシュに変換して精度評価を行う。

本実施例の地上投影変換では、約２０,０００×１０,０００ピクセルの画像を５,０００×５,０００ピクセルの画像へと変換する。観測領域は５,０００×５,０００ｍ²であるため、画像全体の分解能は、１ピクセルあたり１ｍである。内部パラメータは、既知であるとし、表１の値を用いる。ＰＯＣに基づく画像マッチングのウィンドウサイズは、１２８×１２８ピクセルとする。画像マッチングを行う領域は、２つのＳＡＲ画像間で共通する領域でなければならない。そこで、基準点は、５,０００×５,０００ピクセルの画像の共通領域である３,２００×３,２００ピクセルの領域に配置する。また、ＰＯＣ関数のピークによる誤対応除去の閾値を０．１５に設定する。

バンドル調整により最適化を行うパラメータは、航路ｉ（ｉ＝１、２）におけるレーダ入射角θ_i、スラントレンジ方向の分解能の逆数α_vi、外部パラメータＲの回転角度φ、外部パラメータｔのｘおよびｙ成分であるｔ_xおよびｔ_yの７つである。
＜５．３実験結果＞
地上投影変換後のグランドレンジ画像と標高マップを図８に示す。図８（ａ）および（ｂ）は、それぞれ航路１および航路２の地上投影変換後のグランドレンジ画像である。図８（ｃ）は、画像マッチングの際に基準点を配置した領域をグランドレンジ画像から抽出したものである。図８（ｄ）は、バンドル調整前のパラメータを用いて計算した標高マップである。図８（ｅ）は、バンドル調整で最適化されたパラメータを用いて計算した標高マップである。図８（ｆ）は、国土地理院の基盤地図情報（数値標高モデル）データ５ｍメッシュ（標高）を使用した標高マップである。また、本発明の手法によって計算された３次元点群とＤＥＭとの残差と３次元復元点数を表２に示す。

バンドル調整前後のパラメータ値と再投影誤差を表３に示す。以上の結果から、本発明を用いることで１０ｍ以下の誤差で３次元計測が可能であることが確認できる。また、図８（ｄ）、（ｅ）、（ｆ）から、バンドル調整を用いてパラメータの最適化を行うことで、高精度な３次元点群を取得できることがわかる。ＳＡＲ画像の３次元計測においても、コンピュータビジョンで用いられているバンドル調整が有効であることを示している。

本実験で真値として用いたＤＥＭは、人工物や樹木などの影響を除去した地表面の標高モデルである。これに対して、実験で使用したＳＡＲ画像は、樹木などを含めたデータであり、このＳＡＲ画像から求めた３次元点は、樹木の高さに相当する数メートル分の誤差が含まれている。本実験では真値にＤＥＭを用いているが、建物や樹木などの高さを含んだ地表モデルである数値表層モデル（ＤＳＭ：Digital Surface Model）を用いることができれば、提案手法の計測誤差は小さくなると考えられる。

また、計測誤差が生じる原因に、レーダ画像特有のジオメトリック画像変調やスペックルノイズによる画像マッチングへの影響が考えられる。そのため、これらを考慮したＳＡＲ画像のための画像マッチングアルゴリズムを検討する必要がある。たとえば、Dellingerらは、スペックルノイズを考慮したＳＡＲ画像のためのＳＩＦＴアルゴリズムを提案している(非特許文献１０)。

ジオメトリック画像変調の中でも、陰影効果が生じている箇所は、受信信号が０でありデータが存在しない領域であるため、マッチングを正しく行うことができない。表２に示したように最大残差が７７ｍであったが、明らかに誤対応点を用いて３次元計測をした結果であり、全体の計測精度を悪化させる原因となっている。このような影響を除去するためには、画像中の陰影領域を抽出し、陰影領域に基準点を配置しないで画像マッチングを行う必要がある。また、３枚以上のＳＡＲ画像から陰影領域が生じないように適切に画像選択を行うことで、陰影効果による影響を受けずに３次元計測を行えることが容易に予想される。

最も本発明の効果が発揮される利用例として対自然災害の緊急ＳＡＲ観測時における３次元計測が挙げられる。衛星ＳＡＲによるＩｎＳＡＲの場合、周回軌道の制約から３次元計測結果を得るには、通常、数日間以上の日数を要する。例えば陸域観測技術衛星だいちの場合、周回周期の半分の２３日程度以上と考えられる。また、従来のレーダ測量法の場合、高精度な３次元計測結果を得るためには予めＧＣＰ取得の為に災害現場での地上作業が必要とされる上、精密な軌道情報の配信を待つ必要があると考えられる。災害規模によっては現地での地上作業が不可能である事も有り得る事は容易に想像できる。

本発明による３次元計測手法は、精密な軌道情報、及びＧＣＰ取得の為の地上作業を必要としないことから観測画像取得後、即時に３次元計測結果を得る事ができると考えられる。

また、本発明をＩｎＳＡＲと併用することによって、ＩｎＳＡＲの欠点である、位相の２πの不定性を克服することができる。つまり、ＩｎＳＡＲにおける位相の２πの不定性を本発明による結果によって解消することで、正確なＳＡＲ図から正確な立体地形図を形成することができるようになる。

Claims

同一地域のデータを含む少なくとも２枚の互いに異なる航路で取得されたＳＡＲ図を用いて、上記同一地域の所定地点の標高を含む座標を決定する方法であって、
上記所定地点において、地表の対象点の位置のスラントレンジ面への射影における歪みを測定して上記標高を含む座標を決定するものであり、
（１）少なくとも２枚の上記ＳＡＲ図を、同一方角に向き同一縮尺のグランドレンジ面上のＳＡＲ図に変換するステップと、
（２）グランドレンジ面上の上記ＳＡＲ図から共通の領域を選択するステップと、
（３）共通の上記ＳＡＲ図の一方における上記所定地点付近の局所図を動かして、上記ＳＡＲ図の他方との最大一致を示す点を見いだすことで、上記所定地点の異なる上記ＳＡＲ図上での対応点を見いだすステップと、
（４）見いだした各対応点の座標をスラントレンジ面上の座標に変換して、スラントレンジ面上の該対応点の座標の位置から上記所定地点の標高値を見いだすステップと、
を含むことを特徴とする、ＳＡＲ図からの立体地形図形成方法。
上記所定地点のグランドレンジ面上のＳＡＲ図における位置は、上記グランドレンジ面上の上記局所図の重心位置とするものであることを特徴とする請求項１に記載のＳＡＲ図からの立体地形図形成方法。
上記所定地点は、上記同一地域を構成する複数の各点であることを特徴とする請求項１または２のいずれかに記載のＳＡＲ図からの立体地形図形成方法。
同一地域のデータを含む少なくとも２枚の互いに異なる航路で取得されたＳＡＲ図を用いて、上記同一地域の所定地点の標高を含む座標を推定することを含むＳＡＲ図からの立体地形図形成方法であって、
（１）上記ＳＡＲ図の各々にグランドレンジ変換を行って各々のグランドレンジ図を得るステップと、
（２）上記各々のグランドレンジ図について、それらを重ね合わせるために必要な図間の回転角と並行移動量を推定するステップと、
（３）所定の地点に対応する上記各々のグランドレンジ図上の地点の座標を要素とする対応点群を求めるステップと、
（４）上記対応点の各々の要素について、上記グランドレンジ変換の逆変換を行ってスラントレンジ図上の対応点群に座標変換するステップと、
（５）スラントレンジ図上の上記対応点群について、プラットフォーム高度、スラントレンジ方向のビーム入射角、スラントレンジ方向、およびアジマス方向の分解能の逆数を含む上記ＳＡＲ図の内部パラメータと、回転角と並行移動量とを含む外部パラメータと、を用いて上記対応点群の各々の要素の高度を推定するステップと、
（６）上記内部パラメータまたは／および外部パラメータについて、バンドル調整に基づく再投影誤差の最小化を実施し、該再投影誤差が最小となる最適内部パラメータまたは／および最適外部パラメータを求めるステップと、
（７）上記最適内部パラメータまたは／および最適外部パラメータを用いて、高度推定値を得るステップと、
を含むことを特徴とするＳＡＲ図からの立体地形図形成方法。
請求項４の（２）のステップと（３）のステップは、位相限定相関法を用いたステップであることを特徴とするＳＡＲ図からの立体地形図形成方法。
上記ＳＡＲ図は、単偏波スラントレンジ振幅強度図であることを特徴とする、請求項４あるいは５に記載のＳＡＲ図からの立体地形図形成方法。
上記バンドル調整における評価関数は、少なくとも２枚の上記ＳＡＲ図のスラントレンジ面上の対応点群の各要素の座標と上記ＳＡＲの内部パラメータまたは／および外部パラメータを用いて対応点群の各要素の３次元座標を算出し、再度２次元スラントレンジ面に投影した際の対応点群の各要素の座標との差の２乗和であることを特徴とする請求項４から６のいずれか１つに記載のＳＡＲ図からの立体地形図形成方法。