TOP > 国内特許検索 > ランドサットTM画像の精密幾何補正方法及び衛星画像の精密幾何補正方法 > 明細書

明細書 :ランドサットTM画像の精密幾何補正方法及び衛星画像の精密幾何補正方法

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4005336号 (P4005336)
公開番号 特開2003-141507 (P2003-141507A)
登録日 平成19年8月31日(2007.8.31)
発行日 平成19年11月7日(2007.11.7)
公開日 平成15年5月16日(2003.5.16)
発明の名称または考案の名称 ランドサットTM画像の精密幾何補正方法及び衛星画像の精密幾何補正方法
国際特許分類 G06T   1/00        (2006.01)
G06T   3/00        (2006.01)
FI G06T 1/00 285
G06T 3/00 200
請求項の数または発明の数 4
全頁数 16
出願番号 特願2001-342440 (P2001-342440)
出願日 平成13年11月7日(2001.11.7)
新規性喪失の例外の表示 特許法第30条第1項適用 平成13年5月8日 社団法人日本リモートセンシング学会開催の「第30回日本リモートセンシング学会学術講演会」において文書をもって発表
審査請求日 平成16年4月2日(2004.4.2)
特許権者または実用新案権者 【識別番号】503360115
【氏名又は名称】独立行政法人科学技術振興機構
発明者または考案者 【氏名】飯倉 善和
個別代理人の代理人 【識別番号】100091443、【弁理士】、【氏名又は名称】西浦 ▲嗣▼晴
審査官 【審査官】岡本 俊威
参考文献・文献 特開2001-160133(JP,A)
特開平05-187876(JP,A)
木村由貴子 他,幾何補正のためのGCPの自動取得とその評価,第28回学術講演論文集,日本,(社)日本リモートセンシング学会,2000年 5月,p191-194
調査した分野 G06T 1/00
G06T 3/00
特許請求の範囲 【請求項1】
地図上の座標とランドサットTM画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて前記ランドサットTM画像の座標と前記地図の座標とを整合させるために前記ランドサットTM画像を幾何補正するランドサットTM画像の精密幾何補正方法であって、
前記幾何変換式として前記ランドサットTM画像に付随して提供されるシステム情報に含まれる平行移動パラメータを含む幾何変換式を用い、
前記ランドサットTM画像全体に現れる画像上の特徴と前記地図上の地図情報に関連して得られる地図関連情報のうち前記画像上の特徴に対応する対応地図情報との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて、前記平行移動パラメータの最適化を行い、
前記画像上の特徴が、前記画像上に現れる陸地と海との相違であり、
前記対応地図情報が前記地図上に示される陸地と海との相違であり、
前記評価関数として、海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの識別データの正否を、前記対応地図情報により判定して得た識別率を用いることを特徴とするランドサットTM画像の精密幾何補正方法。
【請求項2】
前記システム情報に含まれる幾何変換式は、
p-p=[cosθ*(x-x)-sinθ*(y-y)]/Δ
l-l=[-sinθ*(x-x)-cosθ*(y-y)]/Δ
であり、上記式において、x及びyは前記地図の座標系の中心座標であり、x及びyは前記地図の座標系の任意の位置の座標であり、p及びlは前記ランドサットTM画像の座標系の中心座標であり、p及びlは前記ランドサットTM画像上の任意の位置の座標であり、θは前記地図の座標系に対する前記画像の座標系の回転角であり、Δは前記ランドサットTM画像の空間解像度であり、
上記式を、p=ax-by+c及びl=bx-ay+dと簡略化したときの(但しa乃至dは変数)、前記c及びdを前記平行移動パラメータとすることを特徴とする請求項1に記載のランドサットTM画像の精密幾何補正方法。
【請求項3】
前記システム情報に含まれる幾何変換式は、
p-p=[cosθ*(x-x)-sinθ*(y-y)]/Δ
l-l=[-sinθ*(x-x)-cosθ*(y-y)]/Δ
であり、上記式において、x及びyは前記地図の座標系の中心座標であり、x及びyは前記地図の座標系の任意の位置の座標であり、p及びlは前記ランドサットTM画像の座標系の中心座標であり、p及びlは前記ランドサットTM画像上の任意の位置の座標であり、θは前記地図の座標系に対する前記画像の座標系の回転角であり、Δは前記ランドサットTM画像の空間解像度であり、
前記p及びlが前記評価関数を用いて最適化するパラメータであることを特徴とする請求項1に記載のランドサットTM画像の精密幾何補正方法。
【請求項4】
地図上の座標と衛星画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて前記衛星画像の座標と前記地図の座標とを整合させるために前記衛星画像を幾何補正する衛星画像の精密幾何補正方法であって、
前記衛星画像上に海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が、陸地であるか海であるかの判定の正否を前記地図上に示される陸地と海のデータに基いて判定して、その識別率を求め、
前記識別率を評価関数として前記複数のパラメータの適切値を求め、
前記適切値を前記幾何変換式に入力して前記衛星画像を幾何補正して幾何補正された衛星画像を定め、
数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、
前記シミュレーション画像及び前記幾何補正された衛星画像のそれぞれの画素値の相関係数を評価関数として、前記複数のパラメータの最適値化を行うことを特徴とする衛星画像の精密幾何補正方法。
発明の詳細な説明
【0001】
【発明の属する技術分野】
本発明は、ランドサットTM画像の精密幾何補正方法及び衛星画像の精密幾何補正方法に関するものである。
【0002】
【従来の技術】
ランドサットTMデータ(Thematic Mapper Data)の標準処理データには衛星画像を地図上(UTM座標系)に投影するための情報(システム情報)として、例えば複数のパラメータを含む幾何変換式が含まれている。しかし、システム情報に基づく幾何変換(システム幾何変換)では、かなりの誤差が生じるため、一般にはほとんど用いられることがない。宇宙開発事業団から入手できるランドサットTMデータは、実際の衛星画像に対して所定の幾何補正が施されて変換された画像である。しかしながらこの変換された画像の幾何補正精度は必ずしも高くない。そこで従来、幾何的な精度を要求される場合には、複数のGCP(地上制御点)を利用したアフィン変換に基づく精密幾何補正方法(正射投影を含む)が採用されるのが一般的である。このGCP(地上制御点)を用いた精密幾何補正方法では、目視により同定された複数のGCPの衛星画像上の座標と対応する地図上の座標からアフィン変換式を最小2乗法を用いて統計的に推定する。画素単位のGCPの決定には誤差がさけられないため、精密幾何補正の精度を統計的に上げるには、できるだけ多数のGCPを取得する必要があるとされてきた。しかしGCPの取得は非効率であり、かつ熟練者が行わない場合には高い精度は期待できない問題がある。
【0003】
そこで発明者等は、多数のGCPを自動的に決定する方法として数値標高モデルを用いて作成できる太陽の直達人射照度画像をテンプレート(シミュレーション画像)として、これとシステム幾何変換後の衛星画像を局所的にマッチングさせる方法(シミュレーション画像法)を提案した(「数値標高モデルを用いた衛星画像の地上制御点の同定」と題して「写真測量とリモートセンシング」Vol.37,No.6,1998に発表)。
【0004】
【発明が解決しようとする課題】
発明者が先に提案したシミュレーション画像法を多数のランドサットTM画像に適用した結果、この方法を採用すれば、対象領域がそれほど広くない場合には、システム幾何変換による衛星画像の位置ズレの大きさや方向が、場所にほとんど依存しないことが分かった。これは、恒星センサーを搭載したランドサットの姿勢制御が極めて高い精度で行われているためであると推察される。
【0005】
しかしながらこの従来のシミュレーション画像法では、局所的に衛星画像を地図に投影(マッチング)できる(衛星画像の座標系と地図の座標系を局所的に整合させることができる)ものの、対象領域全体で衛星画像を地図に投影することができなかった。
【0006】
本発明の目的は、多数のGCPを取得することなしに、対象領域全体に対して、精密幾何補正を実行することができるランドサットTM画像の精密幾何補正方法及び衛星画像の精密幾何補正方法を提供することにある。
【0007】
【課題を解決するための手段】
本発明では、システム幾何変換で用いる幾何変換式のパラメータを最適化する。具体的レベルでは、平行移動パラメータを最適化する。このとき衛星画像とシミュレーション画像の相関係数を評価関数として、パラメータを選択(最適化)する。具体的に変化させるパラメータは平行移動パラメータ2個のみであるため、最適解の探索が容易であり、また1画素以内の微調整も可能となる。
【0008】
そこで本発明は、地図上の座標とランドサットTM画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いてランドサットTM画像の座標と地図の座標とを整合させるためにランドサットTM画像を幾何補正するランドサットTM画像の精密幾何補正方法を改良の対象とする。本発明においては、幾何変換式としてランドサットTM画像に付随して提供されるシステム情報に含まれる平行移動パラメータを含む幾何変換式を用いる。
【0009】
例えば、このシステム情報に含まれる幾何変換式は、次の式である。
【0010】
p-p=[cosθ*(x-x)-sinθ*(y-y)]/Δ
l-l=[-sinθ*(x-x)-cosθ*(y-y)]/Δ
上記式において、x及びyは地図の座標系の中心座標であり、x及びyは地図の座標系の任意の位置の座標であり、p及びlはランドサットTM画像の座標系の中心座標(ピクセル番号とライン番号)であり、p及びlはランドサットTM画像上の任意の位置の座標(ピクセル番号とライン番号)であり、θは地図の座標系に対する画像の座標系の回転角であり、ΔはランドサットTM画像の空間解像度である。
【0011】
本発明では、数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、該シミュレーション画像及びランドサットTM画像の画素値の相関係数を評価関数として、適切な平行移動パラメータを求める(評価関数の最適値が得られるパラメータを求める)。ここで太陽の直達入射照度画像とは、太陽の位置情報及び既知の標高情報から得た太陽の直達入射照度の差を含む画像である。太陽の直達入射照度cosβは、cosβ=cosθ*cos(e)+sin(e)*cos(φ-A)の式に基いて演算することができる。但しθは天頂角であり、Aは方位角であり、eが地表面の斜度であり、φが傾斜方位である。画素値とは、画像中の1画素の明るさに比例するデジタル・ナンバーであり、ランドサットTM画像の画素値は、1画素の輝度に比例する値であり、直達入射照度画像の画素値は、1画素の照度に比例する値である。
【0012】
また相関係数は、ランドサットTM画像及びシミュレーション画像中のそれぞれの画素値の関係性を評価する指数であり、-1から+1の間の数値となる。相関係数は、一般的な公式に基いて求めることができる。そして相関係数を評価関数とするためには、まずランドサットTM画像とシミュレーション画像の大まかなズレ量を目視等により求めておき、そのズレ量を中心にして一方の画像を他方の画像に対してピクセル方向及びライン方向に所定量ずつ移動させて、特定の画素領域(例えば画像の中心の100×100の画素の領域)における前述の相関係数を求める。そして、この相関係数が最大になる(最適値になる)移動量を決定して、この移動量から二つの画像のズレ量を求め、このズレ量を補正できるように平行移動パラメータの値を決定する。
【0013】
直接的には、p及びlを評価関数を用いて最適化または適切化するパラメータとすればよい。しかし具体的には、上記二つの式を、下記のように簡略化する。
【0014】
p=ax-by+c
l=bx-ay+d
但しa乃至dは変数である。そしてこの式のc及びdを平行移動パラメータとする。そしてこれらのパラメータを、二つの画像の画素値の相関係数を評価関数として、この評価関数が最適値になるように定め、前述のp及びlを求める。
【0015】
このようにして定めた平行移動パラメータを用いて上記幾何変換式に従い、対象領域のすべての座標について幾何変換を行えば、多数のGCPを取得することなく、対象領域全体にわたって高い精度の幾何補正を行うことができる。
【0016】
なお平行移動パラメータをそれぞれ変化させたときのピクセル方向の相関係数及びライン方向の相関係数の値が最大値となる値を、平行移動パラメータの最適値とすれば、1画素内の微調整が可能になる。
【0017】
本発明の特徴をより抽象的に表現すれば、本発明の特徴は、ランドサットTM画像全体に現れる画像上の特徴(例えば、前述の輝度)と地図上の地図情報に関連して得られる地図関連情報のうち画像上の特徴に対応する対応地図情報(例えば、前述の太陽の直達入射照度)との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて評価を行いながら、平行移動パラメータの最適化を行うことである。ここで「空間的整合性を統計的に偏りなく定量的に評価する評価関数」の一つが、前述のTM画像上の輝度に基く画素値と太陽の直達入射照度画像(シミュレーション画像)上の太陽の直達入射照度に基く画素値との間の相関係数である。このような評価関数としては、他のものを用いることもできる。例えば、TM画像上の特徴が、この画像上に現れる陸地と海との相違であり、対応地図情報が地図上に示される陸地と海との相違であるとした場合には、評価関数として、海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの識別データの正否を、対応地図情報により判定して得た識別率を用いることができる。尚「層別」とは、統計処理において、事前情報(海岸線からの距離)により、母集団をいくつかの部分集団(層またはカテゴリー)に分けることを意味する。識別率を評価関数とする場合にも、この評価関数が最適値または適切値になるようにパラメータを定める。
【0018】
本発明の特徴を、衛星画像に対する精密幾何補正方法として更に抽象的に(一般的に)表現すれば、次のように表現できる。
【0019】
すなわち本発明は、地図上の座標と衛星画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて衛星画像の座標と地図の座標とを整合させるように衛星画像を幾何補正する衛星画像の精密幾何補正方法を対象として、衛星画像全体に現れる画像上の特徴と地図上の地図情報に関連して得られる地図関連情報のうち画像上の特徴に対応する対応地図情報との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて複数のパラメータの最適化を行うことを特徴とするものである。
【0020】
さらにパラメータの決定を容易にするためには、2つの段階を経てパラメータを決定すればよい。まず最初に、ランドサットTM画像(衛星画像)上に海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの判定の正否を地図上に示される陸地と海のデータに基いて判定して、その識別率を求める。そしてこの識別率を評価関数として複数のパラメータの適切値を求める。次にこの適切値を幾何変換式に入力して衛星画像を幾何補正して幾何補正された衛星画像を定める。その後、数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、このシミュレーション画像及び幾何補正された衛星画像のそれぞれの画素値の相関係数を評価関数として、複数のパラメータの最適値化を行う。このようにすると目視によりランドサットTM画像(衛星画像)とシミュレーション画像のズレを確認する必要がなくなるので、パラメータの決定の自動化が容易になる。
【0021】
【発明の実施の形態】
以下図面を参照して本発明の実施の形態を詳細に説明する。図1は、本発明の第1の実施の形態の方法のデータ処理の流れを説明するための図である。この実施の形態では、衛星画像としてランドサットTM画像を扱う。ランドサットTMデータのシステム情報に含まれる標準処理データに付属するリーダファイルの地図投影アンシナリーレコードに関する情報には、撮影シーンのセンター(中心)の基準となる地図上の中心座標即ちUTM座標(x,y)、UTM座標系に対する送信するTM画像の回転角(オリエンテーション角θ)、および空間解像度(Δ)が含まれている。空間解像度は、画像を撮影したセンサが分解できる最も小さい対象物の単位である。これらの情報を用いれば、以下の幾何変換式(1)及び(2)を用いて任意の地図上のUTM座標x,yを衛星画像上の位置(ラインおよびピクセルに)に変換(いわゆる逆変換)することができる。これらの式(1)及び(2)は、システム情報に付属している。
【0022】
p-p=[cosθ*(x-x)-sinθ*(y-y)]/Δ …(1)
l-l0=[-sinθ*(x-x)-cosθ*(y-y)]/Δ …(2)
上記式において、x及びyは地図の座標系の中心座標であり、x及びyは地図の座標系の任意の位置の座標であり、p及びlはランドサットTM画像の座標系の中心座標(ピクセル番号とライン番号)であり、p及びlはランドサットTM画像上の任意の位置の座標(ピクセル番号とライン番号)であり、θはオリエンテーション角であり、Δは空間解像度である。
【0023】
上記式をアフィン変換式で簡略化すれば次のようになる。
【0024】
p=ax-by+c
l=bx-ay+d
但しa乃至dは変数である。この場合c及びdが平行移動パラメータとなる。システム幾何変換では上式の係数はすべてが決定される。本実施の形態の方法では、画像座標上での平行移動を意味するcとdの2つのパラメータを変数と考える。また上記式(1)及び(2)におけるp及びlをパラメータ(変数)と考えても同じである。
【0025】
この実施の形態では、ランドサットTM画像を所定の幾何変換を経て幾何補正したもの(変換された衛星画像)とシミュレーションにより作成した直達入射照度画像との相関係数を幾何補正の評価関数として用いてパラメータの最適値を決定する。相関係数は位置決めが正しい程大きくなるため、精密幾何補正に用いる幾何変換式のパラメータを最適化することが可能となる。最適化は相関係数が大きくなるようにパラメータを探索することによって実現される。
【0026】
なお、この実施の形態では、上記式を用いた幾何変換に正射投影を含めるため、上式で得られるピクセル値に起伏に基づく位置ズレを加えた衛星画像上の見かけの位置のデータを利用している。最適解を求めるために利用する幾何変換には正射投影を組み込むものとする。すなわち、上記(1)及び(2)式に地図座標(x,y)を代入して得られるピクセル値pに起伏に基づく位置ズレΔpを加える。ただし、起伏に基づく位置ズレΔpを厳密に求めるには、地球の曲率を考慮した計算が必要となり、計算機への負担が大きくなる。なお厳密な計算による位置ズレが地球の曲率を考慮しない場合の位置ズレよりも10%程大きくなることを利用することにより計算時間の短縮をはかることができる。また衛星直下点の起伏に基く位置ズレは、原画像における欠損データの軌跡から推定すればよい。なおこの点に関しては、「写真測量とリモートセンシング」のVol.37,No.4の12頁~22頁に掲載された「ランドサットTM画像の正射投影とその評価」と題する論文に掲載されている。
【0027】
シミュレーション画像を得るために必要な太陽の直達入射照度cosβは、太陽位置(天頂角θと方位角A)と地表面の斜度eと傾斜方位φの関数として以下の式から与えられる。
【0028】
cosβ=cosθcos(e)+sinθsin(e)cos(φ-A)
ここで太陽位置(天頂角θと方位角A)のデータについては、標準処理データに付属する太陽位置を利用する。斜度eと傾斜方位φを求めるには数値標高モデルが必要となる。数値標高モデルは、国土地理院から発行されている数値地図50m(標高)を投影変換(UTM座標系)および再配列(30m)して用いる。国土地理院から発行されている数値地図50mメッシュ(標高)は、格子点(メッシュ)が等間隔の緯度と経度により構成されており、サイズが約50mである。したがって、これをランドサットTMの画像処理に利用するには、投影変換(UTM座標系)および再配列(空間解像度30m)が必要になるのである。なおこの点の詳細に関しては、日本リモートセンシング学会誌、Vol.21,No.2の150~157頁に掲載された「数値標高モデルの投影変換に用いる内挿方の評価法」と題する論文に掲載されている。本実施の形態では、この投影変換(UTM座標系)および再配列(空間解像度30m)したものから、斜度eと傾斜方位φを求めている。
【0029】
本実施の形態では、図2に示す岩手県内の5箇所(黒丸の中に1から5の数字を記載した地域)を対象地域とした。各地点毎に切り出す衛星画像の大きさは縦430、横430、空間解像度は30mとして数値標高モデルを作成した。
【0030】
衛星画像は1995年6月12日撮影のランドサットTMデータ(パス107,ロウ32)を用いた。撮影時の太陽高度は58度、方位は北から時計周りの方向に112度であった。この情報を用いて太陽の直達入射照度のシミュレーション画像を作成した。地点1のシミュレーション画像を図3に示す。
【0031】
最初にシステム幾何変換後の衛星画像とシミュレーション画像とを目視で比較することにより、おおまかな位置合わせを行った。その結果ピクセル方向に40画素、ライン方向に27画素の位置ズレが確認された。そこで目視により確認した位置ズレ(ピクセル方向に40画素、ライン方向に27画素)を中心にして、変換後のランドサットTM画像をシミュレーション画像に対して画素単位で(1画素ずつ)ピクセル方向及びライン方向に移動させて、ランドサットTM画像の中心に位置する例えば100×100個の画素の画素値とそれらに対応するシュミレーション画像の画素の画素値との相関係数を求めた。その結果を、下記表1に示す。
【0032】
【表1】
JP0004005336B2_000002t.gif上記表1の横軸の数字はピクセル方向に40画素を中心にして±2画素ずらすことを意味し、縦軸の数字はライン方向に27画素を中心にして±2画素ずらすことを意味する。そして表中の数字は、画素単位で移動したときの画素値の相関係数である。表1の結果が、評価関数を示している。表1からは、ピクセル方向に40画素、ライン方向に27画素ずれた位置における相関係数0.695が最も大きくなることが分かる。すなわち評価関数の最適値は、ピクセル方向に40画素、ライン方向に27画素ずらすことである。この段階で、このようなズレ量を得られるように前述の平行移動パラメータc及びdを決定することができる。
【0033】
図4には地点1の変換後の衛星画像(バンド5)を示した。図3に示したシミュレーション画像と比べて、両者の陰影が一致していることが確認できる。
【0034】
次に、1画素以内の微調整を行うために行う、前述の位置ズレを原点とした幾何変換式からの平行移動パラメータの最適値の決定方法を説明する。平行移動パラメータc,dを変化させて作成した衛星画像(ピクセル方向に40画素、ライン方向に27画素ずらしたランドサットTM画像)の中心領域の100×100の画素の画素値とこれに対応するシミュレーション画像の画素の画素値の相関係数を計算した。その結果の例を図5及び図6に示す。ピクセル方向への移動量を決める平行移動パラメータcを変化させた場合の相関係数の変化を図5に示しており、ライン方向への移動量を決める平行移動パラメータdを変化させた場合の相関係数の変化を図6に示してある。図5及び図6を比較すると、パラメータcを変化させた場合のほうが、パラメータdを変化させた場合よりも変化がより明らかであるが、両者ともきれいな対称性がある。2つのパラメータを2次元的に探索することにより最適なパラメータを1画素よりも一桁小さいオーダーで決定することができる。これは相関係数の計算に非常に多くの点を利用できるため、結果の精度が向上したためであると思われる。すなわちパラメータc及びdを変化させて相関係数が最大値になるときのそれぞれのパラメータc及びdの値(この例では0.2及び0.2)をそれぞれの最適値として定めることができる。なお、相関係数を計算する画像範囲を地上被覆が均一で起伏の比較的大きな範囲とすれば0.8をこえる相関係数が得られる。しかしながら、そのようにしても位置ズレの情報としては大きな変化は見られないため、単純に標本の数を大きくするという意味で相関係数をとる範囲は全域にとるのが好ましい。
【0035】
このようにして最適な平行移動量を図1の各地点1乃至5に対して求めた結果を表2に示す。
【0036】
【表2】
JP0004005336B2_000003t.gif表2において、横軸のc及びdはそれぞれパラメータを示し、rは相関係数の値を示している。選択した5点はかなり広い範囲を占めるが、最適な平行移動パラメータc及びdの違いはピクセル方向、ライン方向ともで1.1画素程度であった。この程度の範囲ではcとdを固定しても1画素以内の、ライン方向で約1画素に満たない。これはシステム情報に基づく幾何補正の精度がかなり良いことを示していると考えられる。また、西側の地点(図1の地点2と5)でピクセル方向のマイナス側へのズレが大きくなり、ライン方向ではプラス側へのズレが大きくなるなど系統的な特徴が見られる。このことは、より大きな範囲を問題にする場合には回転を含めた最適化が必要なことを示唆している。
【0037】
上記の説明では、パラメータc及びdの最適値を求めた。しかし上記式(1)及び(2)中のp及びl即ちシーンのセンターをパラメータとすることができる。この場合には、幾何変換式のシーンのセンター(p、l)を変化(δp、δl)させることにより1画素より小さな単位で位置ズレを調整する。また、正射投影を行うことにより、広い領域で精度の高い出力画像を作成することが期待できる。更に最適解の計算は、IDL(Interactive Data Language)(Ver5.3)に組み込まれているPowellの方法を用いて行うことができる。
【0038】
図1の例では、パラメータの最適値化を行った後に式(1)及び(2)の幾何変換式にパラメータを代入して、ランドサットTM画像の精密幾何変換を行っている。しかしながら図7に示すように、ランドサットTMデータを幾何変換する際に、式(1)及び(2)を用い、最初に適当な平行移動パラメータの値をこれらの式中に挿入して、パラメータの最適値化を実行し、その結果で平行移動パラメータを変更するようにしてもよい。
【0039】
本実施の形態では、ランドサットTMの標準処理データを幾何変換する簡便な方法(修正システム幾何変換法)が提案されている。この修正システム幾何変換法ではシステム情報に含まれる地図投影に関する情報を基本に、平行移動と地上の起伏による位置ズレ(正射投影)を加味した変換を考えている。そして、太陽の直達入射照度のシミュレーション画像と衛星画像の画素値の相関係数を評価関数として平行移動のパラメータを調整する。0.1画素単位の微調整も可能である。従来のGCPを利用した精密幾何補正では、GCPの選択と座標変換式の統計的な推定という2段階の操作が必要となるのに対して、本発明の方法では従来問題とされてきたGCPの選択が省略できる利点がある。
【0040】
上記実施の形態では、2つの平行移動パラメータの最適値化を行ったが、最適値化できるパラメータは平行移動パラメータに限定されるものではなく、他のパラメータについても同様の方法で最適値化を実行できるのは勿論である。
【0041】
図8は、本発明の第2の実施の形態のデータ処理の流れを示す図である。上記の第1の実施の形態では、目視で衛星画像とシミュレーション画像とのズレ量を見てから、パラメータの最適値化を行っている。この実施の形態では、目視によるズレ量の検出に変えて、自動的に大域的なズレ量を検出する。そこでこの実施の形態では、大域的な方法(海と陸の違いを利用した方法)で、大域的なズレ量を検出した後、前述の第1の実施の形態で採用している局所的な方法(太陽の直達入射照度の推定値を利用した方法)でパラメータを決定する。
【0042】
具体的に説明すると、この実施の形態では、最初にランドサットTM画像上の海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの判定の正否を地図上に示される陸地と海のデータに基いて判定して、その識別率を求める。そしてこの識別率を評価関数として平行移動パラメータの適切値を求める。その後適切値を幾何変換式に入力してランドサットTM画像を幾何補正して幾何補正ランドサットTM画像を定める。次に数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、シミュレーション画像及び幾何補正ランドサットTM画像のそれぞれの画素値の相関係数を評価関数として、平行移動パラメータの最適値化を行う。
【0043】
以下に大域的な方法(海と陸の違いを利用した方法)について説明する。衛星画像で識別できるもっとも大きな特徴は、陸域と海域の違いであり、これらは特に近赤外バンドの画像では明瞭に識別できる。地図座標上の各点が海か陸かの判定は国土地理院から発行される数値情報によって知ることができる。この特徴を利用するため、陸域と海域から複数の評価点を選択する。次に設定した幾何補正式を利用して評価点の衛星画像上の座標を求める。この座標上の輝度値から判定できる海と陸が地図情報と一致する評価点の割り合いを適合度とする。この方法では、評価点の設定が重要である。海岸に近いところだけでなく遠いところも含めることによりパラメータの初期値によらない最適化が可能となる。また評価点を調整することにより、計算時間の短縮を計ることができるため、シーン全体などの広い地域を対象にする場合や多数のパラメータを推定する場合に有効である。
【0044】
この実施の形態では、宇宙開発事業団から提供されたランドサット5号からの衛星画像(TMデータ)「1995年6月12日撮影(パス107/ロウ32)」を用いた。ここで使用するデータはいわゆる標準処理データ(レベル1b)といわれるもので、(1)検知素子間の感度差の調整、および(2)衛星の軌道・姿勢情報を用いたUTM座標系への投影と再配列(システム幾何補正)が行われたものである。しかし、システム幾何補正では(1)地表面の起伏の影響が考慮されていない、(2)軌道情報が正確ではない、(3)画像の配列がUTM座標系のxy軸と対応していない(オリエンテーション角)などの問題があるため、地図のデータと正確には重ね合わせることができない。
【0045】
図9に可視バンドの画像ファイルの一例の画像を示している。ファイルの最初にはヘッダレコード(1行)、各レコードの最初の32バイトはレコードヘッダ、最後の68バイトはレコードトレイラであるが、本解析に必要な情報は含まれていない。他の部分が画像データであるが、左右の黒い部分(輝度値0)はシステム補正の際に生じた撮影されていない部分である。精密幾何変換(正射投影)では画像を各ライン毎のこの境界の軌跡の中間点を衛星直下点とする。
【0046】
標準処理データには、画像ファイルの他に情報ファイル(リードファイル)が付属している。このファイルには幾何変換に必要な以下の情報が記載されている。
【0047】
太陽高度 θ= 62 度
太陽方位 φ= 119 度
オリエンテーション角θ= 0.1830542 radian
シーン中心(UTM座標系)
= 534146.1182 m
= 4463628.9062 m
空間分解能 28.5m
また、この画像ファイル上のシーン中心は
= 3492.0 pixel
= 2983.5 line
である。
【0048】
これらの情報からUTM座標(x、y)が与えられた場合の画像上の位置(p、l)は上記第1の実施の形態で用いた(1)式及び(2)式を用いて次の式で求められる。
【0049】
p=p+[cosθ*(x-x)-sinθ*(y-y)]/28.5
l=l+[-sinθ*(x-x)-cosθ*(y-y)]/28.5
図10に地図上に衛星データの範囲を示した。図10において、領域Aが図9の衛星画像の範囲であり、領域Bが大地的な方法を実施する範囲であり、領域Cが局所的な方法を実施する範囲である。
【0050】
図8に示したデータ処理の流れは、(A)評価用のデータの作成と(B)評価関数の最適化による幾何変換式の推定に分けられる。後者には海と陸の判別を用いる大域的な方法と直達入射照度を用いる局所的な方法がある。それぞれ独立に用いることもできるが、ランドサットTMの場合には、大域的な方法で回転を含む幾何変換式を作成した後で、その平行移動パラメータを修正する形で局所的な方法を用いるのが合理的である。
【0051】
図8のデータ処理の流れでは、まず最初に、評価用のデータを作成する。評価用のデータは2種類(海と陸の判別データと太陽の直達入射照度)あるが、どちらも国土地理院発行の数値地図50mメッシュ(標高)から作成する。
【0052】
図11に示した対象とする衛星画像(パス107/ロウ32用)の範囲と海岸線の多くが重なる矩形の領域Bの数値標高モデルを切り出したものが図11に示すものである。この画像を用いて海岸線からの距離画像(海と陸が接するものをレベル1)を作成した。海岸線からの距離が約750mの領域を距離で15のレベルに層別し、各レベル(カテゴリー)毎に200画素(海と陸が半々)を無作為に抽出し、合計3000個の地点を評価点とし表3に示すようなデータを作成した。この一部を示したものが図12である。図12において、+印が評価点を示している。さらに幾何変換式の決定で用いるデータはこの経緯度で示される位置データをUTM座標系に変換する。
【0053】
【表3】
JP0004005336B2_000004t.gif図13は、図10に示した地域Cの太陽の直達入射照度画像(シミュレーション画像)を示している。この画像は、数値標高モデル(UTM座標系/空間分解能30 m)を数値地図50mメッシュ(標高)を投影変換/再配列(共1次内挿法)して作成したものである。画像の大きさは430×430である。太陽の直達入射照度cosβの求め方は第1の実施の形態の場合と同じである。
【0054】
次に評価関数の最適化による幾何変換式の推定について説明する。
【0055】
図9に示したバンド4の衛星画像に対して、大域的な方法を適用する。システム情報に基づいて評価点の衛星画像上の位置を計算し、各評価点の輝度値(画素値)が閾値γ(20を設定)以下の場合を海域とし、輝度値が閾値より大きい場合を陸域と判断した。そしてこの判断結果(識別データ)の正否を、前述の数値地図50mメッシュ(標高)から得た地図関連情報に照らして判断し、識別率を求めた。その結果、評価点から衛星画像の範囲外にある122地点を除いた2878地点の中で718地点が間違って判定された。
【0056】
この例では、間違った数を識別率とし、これを評価関数として評価値が大きくなるように(識別率が小さくなるように)シーンのセンターのズレ(Δp、Δl)および閾値(γ)に関して最適化を行ったところ、間違いは51個に減少した。ここで最適値化は、Δp、Δlおよび閾値γを周知の最適化手法に従って可変し、その都度、輝度値(画素値)と閾値との対比して行った。このときのパラメータの値は以下の通りである。
【0057】
ケースA: Δp=40.90、Δl=27.18、γ=28.39
ケースAは平行移動のみであるが、回転を含めるために、さらに、オリエンテーション角のズレΔθを含めて最適化を行ったところ、間違いは40個となった。図14に間違った場所と範囲外の地点を地図上に示した。北部に見えるのは衛星画像の範囲外となった評価点である。このときのパラメータの値は以下の通りである。
【0058】
ケースB: Δp=40.25、Δl=27.18、Δθ=0.000217、γ=29.08
空間分解能28.5mも変化させて最適化を行った場合には、間違いの現象はみられなかった。なお、間違った地点を個別的にチェックはしていないが、陸奥小川原湖周辺の湖沼の影響が多く見られるほかは、系統的な傾向は見られない。小さな島や河口、あるいは港の回収などの影響によると思われる。
【0059】
前述の大域的な方法で得られた幾何変換式を用いて、画像の再配列を行った。すなわち幾何補正したランドサットTM画像とシミュレーション画像との間の大まかなズレ量を決定した。この画像の再配列が、第1の実施の形態における目視によるシーンのセンターのズレ量の決定に相当する。その後第1の実施の形態と同様にして、画素単位で幾何補正したランドサットTM画像をピクセル方向及びライン方向に移動して、変換された衛星画像とシミュレーション画像の所定領域におけるそれぞれの画素値の相関係数の変化を求めた。表4及び表5は、それぞれ前述のケースA及びBのパラメータを用いた場合に得られた相関係数の変化を示している。なおこの例では、地形の起伏による誤差を避けるために正射投影を考慮した変換を行った。また内挿には共1次内挿法を用いた。解析には変換後の衛星画像をピクセル方向とライン方向にずらして求めた相関係数を用いた。正しい変換であればズレがないところで最大値をとる。
【0060】
【表4】
JP0004005336B2_000005t.gif【表5】
JP0004005336B2_000006t.gif表4および表5には、ともに横軸にピクセル方向への移動量(画素数)を示して縦軸にライン方向への移動量(画素数)を示している。横軸及び縦軸の数字は、ケースA及びケースBで求めたシーンのセンターのズレΔp及びΔlの値を基準にして(0にして)画素単位で画像をピクセル方向及びライン方向にずらすことを意味している。
【0061】
ケースA(平行移動のみ)の場合も、ケースB(回転も含める)の場合も、共にライン方向にはズレがなく、ピクセル方向に最大2画素程度のズレが見られる。両者を比較すると、後者の方の位置ズレが若干小さい。
【0062】
次に相関係数を評価関数として、各ケースについてパラメータの最適化を行った。その結果を以下に示す。
【0063】
ケースA: Δp=38.79、Δl=27.45、相関係数0.692
ケースB: Δp=39.19、Δl=27.21、Δθ=-0.000213、相関係数0.692
最適化を行った場合には、回転を含める効果はほとんど現れない。これは対象画像の範囲(約10km)では回転の影響は数m(10000×Δθ)程度しかあらわれないためである。したがって、回転のズレを求めるには大域的な方法が不可欠である。また、局所的な方法で最適化を行う場合、初期値を適切に設定しないと解がなかなか見つからないことが起きる。大域的な方法で得られたパラメータを初期値として利用することにより、この問題を解決することができる。
【0064】
【発明の効果】
本発明によれば、相関関数を評価関数として、適切な平行移動パラメータを求めることができるので、この平行移動パラメータを用いて対象領域のすべての座標について幾何変換を行えば、多数のGCPを取得することなく、対象領域全体にわたって高い精度の幾何補正を行うことができる利点がある。
【図面の簡単な説明】
【図1】本発明の第1の実施の形態の方法のデータ処理の流れを説明するための図である。
【図2】対象地域を示す図である。
【図3】図2に示した地点1のシミュレーション画像を示す図である。
【図4】図2に示した地点1の変換後の衛星画像(バンド5)を示す図である。
【図5】パラメータcに対する相関係数の変化を示す図である。
【図6】パラメータdに対する相関係数の変化を示す図である。
【図7】第1の実施の形態の別のデータ処理の流れを示す図である。
【図8】本発明の第2の実施の形態のデータ処理の流れを示す図である。
【図9】可視バンドの画像ファイルの画像を示す図である。
【図10】対象地域を示す図である。
【図11】図10の領域Bの数値標高モデルを示す図である。
【図12】無作為に抽出した評価点の例を示す図である。
【図13】図10に示した地域Cの直達入射照度画像(シミュレーション画像)を示す図である。
【図14】陸と海の判定を誤った評価点を示す図である。
【符号の説明】
A乃至C 領域
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13