(57)【特許請求の範囲】
【請求項１】
同一地域の異なる解像度の２種類の画像の組み合わせを用いたイメージングマッチング方法において、
前記２種類の画像に対して標定作業を行って関数で関係づけた拡張ステレオモデルを生成する工程と、
前記拡張ステレオモデルの両方の画像に共通に存在する既知の地物のポリゴン点列を基準の不定形窓とする工程と、
前記不定形窓にそれぞれ高さ（ｚ）を与えて、この高さ（ｚ）を変化させる工程と、
前記不定形窓の高さが変化される毎に、その高さに応じた前記拡張ステレオモデルの両方の画像上における投影位置を求める工程と、
前記拡張ステレオモデルの両画像の投影位置における前記不定形窓で区切られる領域の輝度値を比較する工程と、
前記比較結果が最も近似したときの不定形窓の高さ（ｚ）を前記拡張ステレオモデルの両画像の地物の高さとして出力する工程と、
前記出力された高さ（ｚ）を用いて前記地物の３次元座標を決定する工程と
を有することを特徴とする不定形窓を用いた画像間拡張イメージングマッチング方法。
【請求項２】
地域の地表の地理データを備え、地物の高さが求められる毎に、前記地理データから前記高さを減算して地物の比高を求める工程と
を有することを特徴とする請求項１記載の不定形窓を用いた画像間拡張イメージングマッチング方法。
【請求項３】
前記不定形窓に基づく画像が撮影日が新しい画像に存在しないときは、前記地物が消失していることを知らせる工程と
を有することを特徴とする請求項１又は２記載の不定形窓を用いた画像間拡張イメージングマッチング方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
航空写真や衛星画像など、同一地域を撮像した２種類以上の画像を利用して、標高を持たない２次元的な数値地図データを３次元データ化（標高値を付与）したり、既成の数値地図データに記述された地物の滅失を確認する手法に関する。
【０００２】
【従来の技術】
従来の写真測量では、計測対象を同種のカメラを用いて異なる二箇所の位置（もしくはそれ以上）から撮影し、各画像上で対応点を求め、三角測量の原理に基づいて３次元空間内の位置を決定する作業を行う。
【０００３】
この作業を自動的に行うための方法には、ステレオ画像を構成する一方の画像内の１点もしくは規則的な格子上に区分された領域、あるいは線構造（エッジ）に着目し、それに対応する点もしくは領域、あるいは線構造を、もう一方の画像で探索しマッチングするステレオマッチング方法が一般的である。
【０００４】
【発明が解決しようとする課題】
しかし、この方法では、あらかじめ形状が分かっている領域が存在したとしても、例えば撮影高度、撮影日等が相違しているからその形状を探索条件とすることができない。
【０００５】
そのため、この形状を拘束条件とするマッチングが行えず、効率的な標高の決定ができないことになり、結果としてその地物の３次元データを得ることができない。
【０００６】
これは、例えば既往の数値地図データに記述された地物、例えば建物などに対して標高を与え３次元データ化しようとする際に、ステレオ画像によって計測された点や線を再度、データ編集しなければ３次元数値地図データの生成が行えないと言う問題になる。
【０００７】
また、最近、高分解能商業衛星が実用化されたり、ラインスキャナを航空機に搭載しデジタル画像データを取得する試みも始まり、航空写真一種類に限らず同一地域を異なる撮像装置で画像化したデータが個別に入手される状況が広がり始めている。
【０００８】
すなわち、従来の写真測量のように予めステレオの撮影を意図せず、同一地域を異なる位置から撮影した画像が個別に取得されるようになってきている。
【０００９】
２１世紀に向けた情報化社会においては、地理情報システム（ＧＩＳ）の普及とともに、その情報基盤となる数値地図データをタイムリーに更新してデータを維持管理することが必要となる。
【００１０】
しかし、最新の数値地図をタイムリーに得ようとしても、従来のような数年に１回のステレオ航空写真の撮影によるデータ更新では、適時更新は困難であり、しかもコストも低減できない。
【００１１】
そこで、同一地域を適時、撮像し個別に取得される多種多様な画像を用いて従来の写真測量と同等の計測を行ったり地図更新を行うといったことが有用になってくる。しかし、現実にはそのようなシステムは存在しない。
【００１２】
本発明は、以上の課題を解決するためになされたもので、同一地域を撮影した２種類以上の多種多様な画像を用いて標高を持たない数値地図データを容易に３次元データ化することによって、既存の地図の更新ができる画像間拡張イメージマッチング方法を得ることを目的とする。
【００１３】
【課題を解決するための手段】
本発明は、同一地域の異なる解像度の２種類の画像の組み合わせを用いたイメージングマッチング方法において、前記２種類の画像に対して標定作業を行って関数で関係づけた拡張ステレオモデルを生成する工程と、前記拡張ステレオモデルの両方の画像に共通に存在する既知の地物のポリゴン点列を基準の不定形窓とする工程と、前記不定形窓にそれぞれ高さ（ｚ）を与えて、この高さ（ｚ）を変化させる工程と、前記不定形窓の高さが変化される毎に、その高さに応じた前記拡張ステレオモデルの両方の画像上における投影位置を求める工程と、前記拡張ステレオモデルの両画像の投影位置における前記不定形窓で区切られる領域の輝度値を比較する工程と、前記比較結果が最も近似したときの不定形窓の高さ（ｚ）を前記拡張ステレオモデルの両画像の地物の高さとして出力する工程と、  前記出力された高さ（ｚ）を用いて前記地物の３次元座標を決定する工程とを備えたことを要旨とする。
【００１４】
【発明の実施の形態】
＜実施の形態１＞
本実施の形態１では、図１に示すように、航空機（ヘリコプタを含む）及び衛星で取得した同一地域の航空写真画像、衛星写真画像、ラインスキャナ画像の内の任意の組み合わせ〔衛星画像（写真又はラインスキャナ画像）と航空画像（写真又はラインスキャナ）との組み合わせ又は異なる高度の航空写真同士の組み合わせ若しくは異なる高度のスキャナ画像の組み合わせ等の組み合わせを含む〕を、コンピュータでデジタルフォーマット化してハードディスク等に格納して、この２つの解像度を合わせ（撮影高度、縮尺に基づく）、これらの画像に対して標定作業）を行う。
【００１５】
そして、各々の画像に対して地上の３次元座標を各画像上の２次元座標に投影するための関数を決定する。このように地上の３次元座標を各々の画像座標に投影する関数によって２つの画像が関係づけられた画像の組み合わせを拡張ステレオモデルＤｉ（輝度画像Ｉ１、輝度画像Ｉ２）と呼ぶことにする。
【００１６】
このような拡張ステレオモデルＤｉを用いると、例えば、標高値を持たず２次元的に記述された既往の数値地図上の地物に標高を与え、３次元データ化することができる。
【００１７】
すなわち、標高を求めたい任意の形状をした地物Ｂｕ（既知）を閉領域の基準の不定形窓Ｐｌ（ポリゴン列）とし、これが任意の標高地（高さ方向）をもつように定める。
【００１８】
そして、この不定形窓Ｐｌの高さを変化させながら輝度画像Ｉ１、輝度画像Ｉ２に前述の関数を用いて輝度画像Ｉ１、Ｉ２上の位置を求め、この位置における不定形窓Ｐｌを投影したときの輝度画像Ｉ１、Ｉ２における輝度値を比較し、評価関数によって最も一致（近似）したと判定されたときの不定形的Ｐ１の高さ（標高値Ｈ）をその輝度画像における地物の高さと決定する。図１では標高値がＨ３のときに最も一致しているとした例である。
【００１９】
また、図１では、理解を容易にするために、不定形窓Ｐｌの高さを変化させて輝度画像Ｉ１、Ｉ２に投影したときの不定形窓Ｐｌの形状を変えて示している。
【００２０】
＜実施の形態２＞
また、地図更新において地図にない新築建物などを３次元データとして追加する場合がある。この場合、従来の方法では、オペレータがステレオ観測しながら３次元データ化しなければならないが、異なる種類の画像の組み合わせではオペレータはステレオ観測ができないし、そのようなシステムも存在しない。
【００２１】
しかし、拡張ステレオモデルＤｉを用いると、この拡張ステレオモデルＤｉは輝度画像Ｉ１と輝度画像Ｉ２とが関数によって関係づけられているので、図２に示すように、例えば航空機によって得た輝度画像Ｉ１を画面に表示し、この画像上でマウス操作等で地物Ｂｕの輪郭に相当する閉領域Ｂ１を与える。そして、この閉領域に投影が一致する位置に不定形窓Ｐｌｈを発生させ、標高Ｈ（ｚ）を変化させる。
【００２２】
例えば３０ｍ、８０ｍ、１５０ｍの順に標高を変化させると、３０ｍ、８０ｍ、１５０ｍと変化したときの輝度画像Ｉ１、Ｉ２における不定形窓Ｐｌｈの輝度画像Ｉ１、Ｉ２の投影位置をその都度求めながら、この位置における不定形窓Ｐｌｈを投影したときの輝度画像Ｉ１、Ｉ２における輝度値を比較し、評価関数によって最も一致（近似）したと判定されるときの不定形的Ｐ１ｈの高さ（標高値Ｈ）をその輝度画像における地物の高さと決定する。従って、地物Ｂｕの三次元座標を推定することができている。
【００２３】
また、図２では、理解を容易にするために、不定形窓Ｐｌｈの高さを変化させて輝度画像Ｉ１、Ｉ２に投影したときの不定形窓Ｐｌｈの形状を変えて示している。
【００２４】
つまり、実施の実施の形態１及び実施の形態２では図３に示すように、前述の組み合わせの画像同士を関数で関係させた拡張ステレオモデルの輝度画像Ｉ１、Ｉ２において、共通に写っている地物Ｂｕから得た不定形窓Ｐｌｈの高さを変化させることによって、輝度画像Ｉ１、Ｉ２における投影位置を求め、この投影位置における不定形窓で区切られる輝度値の評価関数の一致（近似）で地物Ｂｕの高さを推定している。
【００２５】
＜実施の形態３＞
また、数値地形モデルが利用できる場合は、図４に示すように、地物位置での標高を地物標高ｚから引くことにより、地物の比高も得る。これによって、正確な建物の高さを得ることができる。
【００２６】
更に、例えば既往の数値地図に記述された建物に着目するとする。この比高が０に近い値に決定されたとすると、建物が滅失したものとして推定することができる。既に３次元座標が付与された建物についても、新規に撮影された画像と旧画像とで地物が投影される位置の画像を比較することにより、地物の滅失や変更を推定することもできる。
【００２７】
次に、上記実施の形態１、２、３の具体例を以下に説明する。図５は、本発明の地物の標高および比高の推定方法の構成図である。図６は図５の標高探索部の概略構成図である。
【００２８】
例えば、デジタル画像入力部１は、２つの画像をデジタル画像の形で画像処理装置であるコンピュータのハードディスク２等に格納する。このとき、航空写真のようにフィルムや印画紙として供給される画像の場合は、スキャナで読み取り、Ａ／Ｄ変換部によってデジタル画像に変換する。この変換した２つの写真のデジタル画像のピクセル座標（ｘ１ｐｉｘ，ｙ１ｐｉｘ）および（ｘ２ｐｉｘ，ｙ２ｐｉｘ）における画像輝度をＩ１（ｘ１ｐｉｘ，ｙ１ｐｉｘ）およびＩ２（ｘ２ｐｉｘ，ｙ２ｐｉｘ）とする。
【００２９】
また、前述のハードディスク２には、拡張ステレオモデル撮影時の対地高度Ｈと、画像間の撮影位置間隔（基線長Ｂ）と、２つの画像の１ピクセルに対応する地上解像度平均Ｄと、既知の建物のポリゴン（不定形窓）等が記憶されている。
【００３０】
次に、標定要素算出部３は、２つの輝度画像Ｉ１、Ｉ２について、地上の３次元座標が投影される各画像の画像座標算出のための座標変換関数ＦＸ、逆座標変換関数関数ＧＸ（総称して標定情報ともいう）を決定する。これは航空写真測量において標定作業と呼ばれているものにあたる。標定作業の方法や求めるパラメータは画像センサ（例えばエリアセンサとリニアアレイセンサ）によって異なるが、最終的には、３次元空間内の点（ｘ，ｙ，ｚ）が投影される左右画像のピクセル座標（ｘ１ｐｉｘ，ｙ１ｐｉｘ）および（ｘ２ｐｉｘ，ｙ２ｐｉｘ）を計算する関数を決定できる。
【００３１】
つまり、次の４つの座標変換関数ＦＸ１（ｘ，ｙ，ｚ）、ＦＹ１（ｘ，ｙ，ｚ）、ＦＸ２（ｘ，ｙ，ｚ）、ＦＹ２（ｘ，ｙ，ｚ）を解析的、もしくは手続き的に決定できる。
【００３２】
ｘ１ｐｉｘ＝ＦＸ１（ｘ，ｙ，ｚ）            式（１）
ｙ１ｐｉｘ＝ＦＹ１（ｘ，ｙ，ｚ）            式（２）
ｘ２ｐｉｘ＝ＦＸ２（ｘ，ｙ，ｚ）            式（３）
ｙ２ｐｉｘ＝ＦＹ２（ｘ，ｙ，ｚ）            式（４）
また、ｚ＝一定として式（１），式（２）を解析的又は手続的に解くと、各画像上の点がある標高ｚにあると仮定した時の地上でのｘ座標とｙ座標を求めることができる。
【００３３】
ｘ＝ＧＸ１（ｘ１ｐｉｘ，ｙ１ｐｉｘ，ｚ）          式（５）
ｙ＝ＧＹ１（ｘ１ｐｉｘ，ｙ１ｐｉｘ，ｚ）          式（６）
同様に式（３），式（４）について解くと、
ｘ＝ＧＸ２（ｘ２ｐｉｘ，ｙ２ｐｉｘ，ｚ）          式（７）
ｙ＝ＧＹ２（ｘ２ｐｉｘ，ｙ２ｐｉｘ，ｚ）          式（８）
である。
【００３４】
航空写真やエリアＣＣＤセンサの様なエリアセンサ画像の場合については、関数ＦＸ１（ｘ，ｙ，ｚ）、ＦＹ１（ｘ，ｙ，ｚ）、ＦＸ２（ｘ，ｙ，ｚ）、ＦＹ２（ｘ，ｙ，ｚ）の詳細を図７に記述する。
【００３５】
なお標定作業の方法については既存の方法を用いるものとする。例えばエリアＣＣＤセンサの場合は文献［１］、リニア型中心投影センサについては文献［２］に詳しく記述されている。
【００３６】
［１］日本写真測量学会：解析写真測量改訂版，１９８９。「２」内田修：ステレオ衛星画像を用いた３次元計測の自動化に関する研究，１９８９。
【００３７】
また、図７における例えばエリアセンサ画像の座標変換関数は、
中心投影画像の投影中心の座標をＰ０＝ｔ［Ｘ０  Ｙ０  Ｚ０］、カメラの地上座標に対する回転角をＴ＝ｔ［ω  φ  κ］、これらを合わせたもの（外部標定要素）をＥ＝ｔ［ｔＴ  ｔＰ０］とする（左肩のｔは転置行列もしくは転置ベクトルを示す）。このとき、地上座標Ｐ＝ｔ［ｘ  ｙ  ｚ］はカメラ座標系（投影中心に相対的な３次元座標系）ｐｒで次のように表現される。
【００３８】
ｐｒ＝ｔ［ｘｒ  ｙｒ  ｚｒ］＝Ｒ（Ｔ）（Ｒ−Ｐ０）          式（２２）
ここで、Ｒ（Ｔ）は、次の式で定義される回転行列である。
【００３９】
【数１】
0003618649.tif 000002
また、Ｐが投影される写真座標（投影中心の像（主点）を原点とする２次元座標）ｐ（Ｐ，Ｅ）＝ｔ［ｘ  ｙ］は次の式で表わされる。
【００４０】
【数２】
0003618649.tif 000003
ここでｃは航空写真の焦点距離である。式（２４）は、カメラの外部標定要素と地上点の３次元座標がわかれば、それに対応する写真座標がわかることを意味する。写真座標ｐは、物理的にフィルム面やＣＣＤ素子の上の実寸の座標であるが、フィルムを読みとる際のスキャナのひずみや非線形やＣＣＤ素子の配置誤差が無視できるとすれば、画像座標ｐｉｍｇとの関係は２次元の射影変換Ｈｉｎで表すことができる。以上より、ＥおよびＨｉｎを知ることができれば、地上点の３次元座標から画像座標に変換することが可能となり、４つの座標変換関数ＦＸ１（ｘ，ｙ，ｚ）、ＦＹ１（ｘ，ｙ，ｚ）、ＦＸ２（ｘ，ｙ，ｚ）、ＦＹ２（ｘ，ｙ，ｚ）が求められたこととなる。
【００４１】
また、写真座標ｐ＝（ｘ，ｙ）であるとすると、これはカメラ座標ではｐｒｐ＝（ｘ，ｙ，−ｃ）である。これにより、ｐ＝（ｘ，ｙ）に写っている点の標高Ｚｐがあるとすると、この点の地上座標系における３次元座標Ｐｐ＝（Ｘ，Ｙ，Ｚｐ）は次の式で与えられる。
【００４２】
【数３】
0003618649.tif 000004
画像座標ｐｉｍｇに逆射影変換Ｈｉｎ−１を行った後、式（２５）を適用すれば、Ｚ＝Ｚｐの拘束の元で画像座標から３次元座標への変換を行うことができる。つまり、座標変換関数ＧＸ１（ｘ，ｙ，ｚ）、ＧＹ１（ｘ，ｙ，ｚ）、ＧＸ２（ｘ，ｙ，ｚ）、ＧＹ２（ｘ，ｙ，ｚ）を求めることができる。
【００４３】
一方、数値地図入力部４は、既存の数値地図を入力する。この数値地図は地物（たとえば建物データ）をベクター形式（座標点列データ）で保持しているものとする。数値地図が画像データの場合は、ラスター／ベクター変換によりベクター形式に直す。
【００４４】
次に、処理対象地物データ抽出部５は、オペレータによる指示に従って数値地図データより、標高付けの対象となる地物データを抽出する。数値地図では一般に対象の種類（建物・道路等）の違いがコード付けされているので、抽出の際は対象となるコードの図形のみを抽出すればよい。
【００４５】
そして、抽出した地物Ｂｕをポリゴン点列Ｐｌｉ（閉多角形、一般に最初と最後が同じ点となる点列（ｐ１，ｐ２，…，ｐｎ，ｐ１）で表される）に変換する。  このとき、対象とする地物Ｂｕが建物データの場合は普通ポリゴン点列になっているものが多いので問題はない。ポリゴン点列になっていない地物の場合は、ＣＡＤソフト等で編集してあらかじめポリゴン点列化しておく。
【００４６】
また、画像上地物計測部６は、拡張ステレオモデルＤｉの画像上に建物等の地物Ｂｕがあるが、数値地図データに対象とする地物形状が記述されていない場合には、拡張ステレオ画像Ｄｉ（輝度画像Ｉ１、Ｉ２）のどちらかを画面に表示し、その上でマウス操作により地物形状を入力する。この場合、地物形状はポリゴン点列Ｐｌｉとし、画像座標を列挙した閉図形（不定形窓）として入力する。例えば、図２においては、輝度画像Ｉ１からＢ１を選択し、このＢ１の不定形窓を得る。
【００４７】
次に、標高探索部７は、標高の探索では、画像上地物計測部６で得られた各地物ポリゴン点列Ｐｌｉ（ｉ＝１，…，ｎｉ）について、以下の処理を繰り返す。
【００４８】
標高探索部７における標高探索範囲および探索間隔設定部７ａは、標高の探索範囲（Ｈｍｉｎ，Ｈｍａｘ）および探索間隔ｄＨを設定する。
【００４９】
このＨｍｉｎには、その地域の地表面の標高もしくはそれより小さい値を、Ｈｍａｘにはその地図の建物階数等から予想される地物の最大標高もしくはそれ以上の値を与える。
【００５０】
地形モデル（国土地理院発行の数値地図５０ｍメッシュ（標高）など）が利用できる場合は各地物位置での地表面の標高を推定することができる。
【００５１】
また、利用できない場合は、すべての地物について適当に同じＨｍｉｎとＨｍａｘを設定すればよい。またｄｚには、拡張ステレオモデルから測定される高さ方向の解像度を与えればよい。
【００５２】
たとえば、各画像間の撮影位置間隔（基線長）をＢ、拡張ステレオモデル撮影時の対地高度の平均をＨ、２つの画像の１ピクセルに対応する地上解像度平均をＤとすれば、ｄＨは次の式で与えられる。
【００５３】
【数４】
0003618649.tif 000005
またこの時、探索標高群Ｈｋを次のように定義する。
【００５４】
Ｈｋ＝Ｈｍｉｎ＋ｄＨ×ｋ（ｋ＝０，１，２，…ｎｋ）    式（１０）
ただし、ｎｋはＨｋ≦Ｈｍａｘとなる最大のｋである。
【００５５】
そして、不定形窓定義部７ｂは、地物データのそれぞれのポリゴン点列の内部をステレオマッチングの際の不定形窓として定義する。具体的には、図８に示すように、地物ポリゴン点列内に等間隔のメッシュ（ｄｘ、ｄｙ）を発生させ、これを不定形窓Ｐｌｉとする。
【００５６】
この処理は、ポリゴン点列Ｐｌｉが数値地図（座標が地上座標）から得られたか、それとも画像上の計測結果（座標が画像座標）から得られたかによって処理が異なる。
【００５７】
・数値地図から得られたポリゴン点列の時
地物データが数値地図起源の時は、図８に示すように、地物データＢｕのそれぞれのポリゴン点列の周辺に等間隔のメッシュ点群ＭＰｉを発生させる。メッシュは任意の地上座標に原点（ｏｘ，ｏｙ）を設定し、ｘ方向およびｙ方向にｄｘおよびｄｙの間隔で発生させ、Ｍ（ｏｘ，ｏｙ，ｄｘ，ｄｙ）とする。
【００５８】
なお、ｏｘおよびｏｙは任意に取ってかまわない。またｄｘおよびｄｙは、デジタル航空写真の１ピクセル〜数ピクセルの地上解像度に設定する。
【００５９】
ここでｉ番目の地物Ｆｉを表すポリゴン点列Ｐｌｉで囲まれる領域Ｐｌｉｎｉ内部に存在するメッシュの交点群ＭＰｉを次のように定義する。
【００６０】
ＭＰｉ＝Ｐｌｉｎｉ∩Ｍ（ｏｘ，ｏｙ，ｄｘ，ｄｙ）          式（１１）
このようにして得た不定形窓Ｐｌｉを用いて、例えば図１では、不定形窓Ｐｌの高さを変化させ輝度画像Ｉ１、輝度画像Ｉ２の投影位置を求め、この投影位置における両方の画像の輝度値を比較する。
【００６１】
そして、輝度値の評価関数が最も一致（近似）したときの標高値Ｈを地物Ｂ１、Ｂ２の標高値ｚと決定する。
【００６２】
・画像上の計測結果から得られたポリゴン点列の時は、不定形窓定義部７ｃ（第２の不定形窓定義部ともいう）は、Ｐｌｉが画像上の計測結果の時は、評価関数計算ループ部７ｄの中できめられる標高Ｈｋにあわせてメッシュを作成する。
【００６３】
今、Ｐｌｉが画像１上で計測されたと仮定し、ｉ番目の画像座標を（ｐｘｊ，ｐｙｊ）とすると、この点が標高Ｈｋにあったときの地上
座標（Ｘｊ，Ｙｊ）は、式（５）および式（６）を用いて、
Ｘｊ＝ＧＸ１（ｐｘｊ，ｐｙｊ，Ｈｋ）        式（１２）
Ｙｊ＝ＧＹ１（ｐｘｊ，ｐｙｊ，Ｈｋ）        式（１３）
で表される。こうしてすべてのポリゴン点列上の点を地上座標に戻した後、数値地図の場合と同様に変換して、メッシュ交点群ＭＰｉｋを得る。Ｐｌｉが画像上で計測された場合も式（１２）および式（１３）によって同様にメッシュ点群ＭＰｉｋを得る。
【００６４】
このような不定形窓を用いることによって、例えば図２に示すように、航空機によって得た輝度画像Ｉ１を画面に表示し、この画像上でマウス操作等で座標が既知の地物Ｂ１の輪郭に相当する閉領域（不定形窓Ｐｌｈ：ポリゴン）を与える。そして、この不定形窓Ｐｌｈを地物Ｂｕにおいて、標高Ｈ（ｚ）を変化させる。
【００６５】
前述の評価関数計算ループ部７ｄは、すべてのＨｋ（ｋ＝０，１，２，…，ｎｋ）について、評価値Ｅ（Ｈｋ）を計算する。
【００６６】
【数５】
0003618649.tif 000006
ここで、Ｅｋｊは、メッシュ交点群ＭＰｉのｊ番目の点（ｘｊ，ｙｊ）が標高Ｈｋの高さにあると仮定した時の評価関数値である。評価値Ｅ（Ｈｋ）は、拡張ステレオモデル上の画像の画素値の一致の程度を評価し、一致が高くなる時に小さくなるものとする。
【００６７】
すなわち、標高Ｈｋの高さにあるメッシュ交点群ＭＰｉのｊ番目の点（ｘｊ，ｙｊ，Ｈｋ）が二つの画像Ｉ１およびＩ２の上へ、各々式（１）〜式（４）により投影されるとき、その位置にある画素値Ｖ１ｋｊおよびＶ２ｋｊを
0003618649.tif 000007
ここで、Ｖ１ｋｊ（ｘｊ，ｙｊ）、Ｖ２ｋｊ（ｘｊ，ｙｊ）は投影したときの画素値であり、（ｘｊ，ｙｊ，Ｈｋ）はＨにあるときの不定形窓のポリゴン点群を示す。
【００６８】
上記（１５）式のように表せば、画素値の差の自乗Ｅｄｉｆの総和や、画素値の相関係数に（−１）をかけた値Ｅｃｏｒなどを用いることができる。
【００６９】
例えば、Ｅｄｉｆｋは次の式で表される。
【００７０】
【数６】
0003618649.tif 000008
また、Ｅｃｏｒｋは次の式で表される。
【００７１】
【数７】
0003618649.tif 000009
そして、地物標高決定部７ｅは、ｉ番目の地物Ｆｉの標高ＨＦｉを、評価関数Ｅ（Ｈｋ）を最小にする高さＨｋｍｉｎとする。すなわち、
0003618649.tif 000010
次に、地物高さ計算部７ｆは、地形モデルがある場合は、２次元のポリゴン点列Ｐｌｉの位置の地表面高さを標高値ＨＦｉから引くことにより、地物比高ＨＲｉとして出力する。また、後で信頼性の解析に用いるため、標高値ＨＦｉでの評価関数値Ｅｉも地物の属性情報として出力する。
【００７２】
そして、３次元地物データ出力部７ｅは、２次元のポリゴン点列Ｐｌｉのｚ座標に標高値ＨＦｉを付け加え、更に地表面比高ＨＲｉを属性データとして付け加えた３次元地物データＰｌ３Ｄｉとして出力する。
【００７３】
次に、３次元地図出力部８は、すべてのＰｌ３Ｄｉをあわせて、標高付きの３次元地図として出力する。
【００７４】
次に、標高推定の信頼性の評価について説明を補充する。
【００７５】
次のような場合は、標高推定の信頼性の推定方法として、次の様な例が考えられる。
【００７６】
・評価関数値Ｅｉが大きい時。つまり、あるしきい値をＴ２として、
Ｅｉ≧Ｔ２                            式（１９）
Ｔ２の値の決定方法としては、評価関数値Ｅｉの平均値より、１〜２σ大きい値を設定することが合理的である。
【００７７】
・求めた高さＨＦｉが探索範囲の両端になる場合。つまり、
ＨＦｉ＝Ｈｍｉｎ  ｏｒ  ＨＦｉ＝Ｈｍａｘ          式（２０）
この場合は、探索範囲が正しく設定されなかった可能性がある。
【００７８】
また、地物の滅失や変更の推定する場合は、次のような場合、地物の滅失や変更の可能性が高いと推定できる。
【００７９】
・地表面比高が小さい時。つまり、０に近いしきい値をＴ１として、
｜ＨＲｉ｜≦Ｔ１                      式（２１）
また、信頼性が低くなる原因として、地物の高さ、大きさ等に変更が起こっている場合がある。
【００８０】
例えば、図９に示すように、Ａ年度の撮影時においては地物の標高がＨｏｒｇであっても、数年後（経年変化）に撮影した場合は、標高がＨｏｒｇからＨｎｅｗに変化している場合がある。
【００８１】
また、建物の大きさ（平面形状）も時間の経過に伴って変化している場合もある（Ｃ年度撮影）。
【００８２】
標高の既知の地物の異動の推定の場合は、すでに標高が既知の地物の異動は次のように推定する。
【００８３】
地物の標高をＨｏｒｇとし、標高の探索範囲を
Ｈｏｒｇ−ｄＨ≦Ｈｋ≦Ｈｏｒｇ＋ｄＨにとる。
【００８４】
経年変化がある場合は、推定した標高がＨｏｒｇと大きく異なるか、評価関数値から標高推定の信頼性が低くなる。
【００８５】
また、旧・新の撮影時期の異なる画像を拡張ステレオモデルに構成し、標高推定を行い、標高推定の信頼性が低くなる場合に経年変化があると判定することも可能である。
【００８６】
また、例えば図１０に示すように、地物が複雑な形状であったり、実際には一定の高さではない場合にも対応可能である。建物について実際に航空写真測量やレーザスキャナで計測した結果と比べてみると、建物の最頻値標高に近い値を示すことが実験で確認されている。
【００８７】
また地物高さが０に近いときは、自動的に地物の滅失箇所候補として抽出し、修正の効率化を図ることができる。
【００８８】
＜実施の形態４＞
次に地図の更新について図１１のフローチャートを用いて説明する。本説明では予め３次元地図（旧）と、２次元地図（旧）と、新規撮影の拡張ステレオモデルとがディスクに記憶されているとする。
【００８９】
このような状態において、二次元地図（旧）を読み込むと共に、拡張ステレオモデルとを読み込み、上記説明のように、既知の地物の不定形窓に高さを与え、これを変化させて拡張ステレオモデルの両方の画像上の投影位置を求め、この投影位置における不定形窓で区切られた領域の輝度値をそれぞれ求めて輝度値を比較し、比較結果が最も一致するときの不定形窓の高さを拡張ステレオモデルにおける地物の高さとする拡張ステレオマッチングによる標高付与処理を行い（Ｓ１）、この３次元地図データを途中成果として記憶する（Ｓ２）。
【００９０】
そして、３次元地図（旧）と、この途中成果の３次元地図データと新規の拡張ステレオモデルとから変化を求める拡張ステレオマッチングによる変化抽出処理を行い（Ｓ３）、これを写真上（拡張ステレオモデル）に割り付ける前の途中成果の３次元地図として記憶する（Ｓ４）。
【００９１】
次に、写真上での編集処理を行う（Ｓ５）。これは、デジタル写真測量ワークステーションを用いて行うのが望ましい。
【００９２】
この編集処理は、ステップＳ４で得た途中成果の地物の３次元データを新規の拡張ステレオモデルの上にオーバレイしてモニターに表示し、３次元データと新規に撮影した拡張ステレオモデルとを比較することによってオペレータがマニュアルで編集を行う。ここでは拡張ステレオモデルとしては新規に撮影した一般的なステレオ画像を用いてもよいが、旧２次元地図を作成したときの画像と新規に撮影した画像とから構成されるステレオモデルを使用してもよい。
【００９３】
そして、この編集処理で得られた３次元データに旧の３次元データを更新する処理を行う（Ｓ６）。
【００９４】
つまり、この一連の処理によって、標高を持たない数値地図データを容易に３次元データ化がなされ、かつ既存の地図の更新ができることになる。
【００９５】
【発明の効果】
以上のように本発明によれば、２次元地図を効率的に３次元地図化することができる。特に航空写真を新しくとる必要がない場合は、コンピュータグラフィックスや都市３次元シミュレーションなどの用途に用いての都市３次元モデルを既存地図から廉価に作成することが可能となる。
【００９６】
また、数値地図を３次元化することにより、２次元地図の経年変化修正業務において、新規の撮影した地図をステレオ航空写真上に重ねることができ、修正が容易になる。
【図面の簡単な説明】
【図１】本実施の形態１の不定形窓を用いた拡張イメージング方法の概略構成図である。
【図２】実施の形態２の不定形窓を用いた拡張イメージング方法の概略構成図である。
【図３】不定形窓の高さ方向の変更を説明する説明図である。
【図４】実施の形態３の比高を説明する説明図である。
【図５】各実施の形態の具体的な構成図である。
【図６】標高探索部の構成図である。
【図７】エリアセンサの幾何学的特性を説明する説明図である。
【図８】地物データのメッシュ化の説明図である。
【図９】撮影時期の異なるときの地物データの変化を説明する説明図である。
【図１０】複雑な形状の建物の標高の求め方を説明する説明図である。
【図１１】実施の形態３の地図の更新を説明するフローチャートである。
【符号の説明】
１  デジタル画像入力部
２  ハードディスク
３  標定要素算出部
４  数値地図入力部
５  処理対象地物データ抽出部
６  画像上地物計測部
７  標高探索部
【図１】
0003618649.tif 000011
【図２】
0003618649.tif 000012
【図３】
0003618649.tif 000013
【図４】
0003618649.tif 000014
【図５】
0003618649.tif 000015
【図６】
0003618649.tif 000016
【図７】
0003618649.tif 000017
【図８】
0003618649.tif 000018
【図９】
0003618649.tif 000019
【図１０】
0003618649.tif 000020
【図１１】
0003618649.tif 000021