(57)【特許請求の範囲】
【請求項１】
地図上の座標とランドサットＴＭ画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて前記ランドサットＴＭ画像の座標と前記地図の座標とを整合させるために前記ランドサットＴＭ画像を幾何補正するランドサットＴＭ画像の精密幾何補正方法であって、
前記幾何変換式として前記ランドサットＴＭ画像に付随して提供されるシステム情報に含まれる平行移動パラメータを含む幾何変換式を用い、
前記ランドサットＴＭ画像全体に現れる画像上の特徴と前記地図上の地図情報に関連して得られる地図関連情報のうち前記画像上の特徴に対応する対応地図情報との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて、前記平行移動パラメータの最適化を行い、
前記画像上の特徴が、前記画像上に現れる陸地と海との相違であり、
前記対応地図情報が前記地図上に示される陸地と海との相違であり、
前記評価関数として、海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの識別データの正否を、前記対応地図情報により判定して得た識別率を用いることを特徴とするランドサットＴＭ画像の精密幾何補正方法。
【請求項２】
前記システム情報に含まれる幾何変換式は、
ｐ−ｐ０＝［ｃｏｓθ０＊（ｘ−ｘ０）−ｓｉｎθ０＊（ｙ−ｙ０）］／Δ
ｌ−ｌ０＝［−ｓｉｎθ０＊（ｘ−ｘ０）−ｃｏｓθ０＊（ｙ−ｙ０）］／Δ
であり、上記式において、ｘ０及びｙ０は前記地図の座標系の中心座標であり、ｘ及びｙは前記地図の座標系の任意の位置の座標であり、ｐ０及びｌ０は前記ランドサットＴＭ画像の座標系の中心座標であり、ｐ及びｌは前記ランドサットＴＭ画像上の任意の位置の座標であり、θ０は前記地図の座標系に対する前記画像の座標系の回転角であり、Δは前記ランドサットＴＭ画像の空間解像度であり、
上記式を、ｐ＝ａｘ−ｂｙ＋ｃ及びｌ＝ｂｘ−ａｙ＋ｄと簡略化したときの（但しａ乃至ｄは変数）、前記ｃ及びｄを前記平行移動パラメータとすることを特徴とする請求項１に記載のランドサットＴＭ画像の精密幾何補正方法。
【請求項３】
前記システム情報に含まれる幾何変換式は、
ｐ−ｐ０＝［ｃｏｓθ０＊（ｘ−ｘ０）−ｓｉｎθ０＊（ｙ−ｙ０）］／Δ
ｌ−ｌ０＝［−ｓｉｎθ０＊（ｘ−ｘ０）−ｃｏｓθ０＊（ｙ−ｙ０）］／Δ
であり、上記式において、ｘ０及びｙ０は前記地図の座標系の中心座標であり、ｘ及びｙは前記地図の座標系の任意の位置の座標であり、ｐ０及びｌ０は前記ランドサットＴＭ画像の座標系の中心座標であり、ｐ及びｌは前記ランドサットＴＭ画像上の任意の位置の座標であり、θ０は前記地図の座標系に対する前記画像の座標系の回転角であり、Δは前記ランドサットＴＭ画像の空間解像度であり、
前記ｐ０及びｌ０が前記評価関数を用いて最適化するパラメータであることを特徴とする請求項１に記載のランドサットＴＭ画像の精密幾何補正方法。
【請求項４】
地図上の座標と衛星画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて前記衛星画像の座標と前記地図の座標とを整合させるために前記衛星画像を幾何補正する衛星画像の精密幾何補正方法であって、
前記衛星画像上に海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が、陸地であるか海であるかの判定の正否を前記地図上に示される陸地と海のデータに基いて判定して、その識別率を求め、
前記識別率を評価関数として前記複数のパラメータの適切値を求め、
前記適切値を前記幾何変換式に入力して前記衛星画像を幾何補正して幾何補正された衛星画像を定め、
数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、
前記シミュレーション画像及び前記幾何補正された衛星画像のそれぞれの画素値の相関係数を評価関数として、前記複数のパラメータの最適値化を行うことを特徴とする衛星画像の精密幾何補正方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、ランドサットＴＭ画像の精密幾何補正方法及び衛星画像の精密幾何補正方法に関するものである。
【０００２】
【従来の技術】
ランドサットＴＭデータ（Thematic Mapper Data)の標準処理データには衛星画像を地図上（ＵＴＭ座標系）に投影するための情報（システム情報）として、例えば複数のパラメータを含む幾何変換式が含まれている。しかし、システム情報に基づく幾何変換（システム幾何変換）では、かなりの誤差が生じるため、一般にはほとんど用いられることがない。宇宙開発事業団から入手できるランドサットＴＭデータは、実際の衛星画像に対して所定の幾何補正が施されて変換された画像である。しかしながらこの変換された画像の幾何補正精度は必ずしも高くない。そこで従来、幾何的な精度を要求される場合には、複数のＧＣＰ（地上制御点）を利用したアフィン変換に基づく精密幾何補正方法（正射投影を含む）が採用されるのが一般的である。このＧＣＰ（地上制御点）を用いた精密幾何補正方法では、目視により同定された複数のＧＣＰの衛星画像上の座標と対応する地図上の座標からアフィン変換式を最小２乗法を用いて統計的に推定する。画素単位のＧＣＰの決定には誤差がさけられないため、精密幾何補正の精度を統計的に上げるには、できるだけ多数のＧＣＰを取得する必要があるとされてきた。しかしＧＣＰの取得は非効率であり、かつ熟練者が行わない場合には高い精度は期待できない問題がある。
【０００３】
そこで発明者等は、多数のＧＣＰを自動的に決定する方法として数値標高モデルを用いて作成できる太陽の直達人射照度画像をテンプレート（シミュレーション画像）として、これとシステム幾何変換後の衛星画像を局所的にマッチングさせる方法（シミュレーション画像法）を提案した（「数値標高モデルを用いた衛星画像の地上制御点の同定」と題して「写真測量とリモートセンシング」Ｖｏｌ．３７，Ｎｏ．６，１９９８に発表）。
【０００４】
【発明が解決しようとする課題】
発明者が先に提案したシミュレーション画像法を多数のランドサットＴＭ画像に適用した結果、この方法を採用すれば、対象領域がそれほど広くない場合には、システム幾何変換による衛星画像の位置ズレの大きさや方向が、場所にほとんど依存しないことが分かった。これは、恒星センサーを搭載したランドサットの姿勢制御が極めて高い精度で行われているためであると推察される。
【０００５】
しかしながらこの従来のシミュレーション画像法では、局所的に衛星画像を地図に投影（マッチング）できる（衛星画像の座標系と地図の座標系を局所的に整合させることができる）ものの、対象領域全体で衛星画像を地図に投影することができなかった。
【０００６】
本発明の目的は、多数のＧＣＰを取得することなしに、対象領域全体に対して、精密幾何補正を実行することができるランドサットＴＭ画像の精密幾何補正方法及び衛星画像の精密幾何補正方法を提供することにある。
【０００７】
【課題を解決するための手段】
本発明では、システム幾何変換で用いる幾何変換式のパラメータを最適化する。具体的レベルでは、平行移動パラメータを最適化する。このとき衛星画像とシミュレーション画像の相関係数を評価関数として、パラメータを選択（最適化）する。具体的に変化させるパラメータは平行移動パラメータ２個のみであるため、最適解の探索が容易であり、また１画素以内の微調整も可能となる。
【０００８】
そこで本発明は、地図上の座標とランドサットＴＭ画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いてランドサットＴＭ画像の座標と地図の座標とを整合させるためにランドサットＴＭ画像を幾何補正するランドサットＴＭ画像の精密幾何補正方法を改良の対象とする。本発明においては、幾何変換式としてランドサットＴＭ画像に付随して提供されるシステム情報に含まれる平行移動パラメータを含む幾何変換式を用いる。
【０００９】
例えば、このシステム情報に含まれる幾何変換式は、次の式である。
【００１０】
ｐ−ｐ０＝［ｃｏｓθ０＊（ｘ−ｘ０）−ｓｉｎθ０＊（ｙ−ｙ０）］／Δ
ｌ−ｌ０＝［−ｓｉｎθ０＊（ｘ−ｘ０）−ｃｏｓθ０＊（ｙ−ｙ０）］／Δ
上記式において、ｘ０及びｙ０は地図の座標系の中心座標であり、ｘ及びｙは地図の座標系の任意の位置の座標であり、ｐ０及びｌ０はランドサットＴＭ画像の座標系の中心座標（ピクセル番号とライン番号）であり、ｐ及びｌはランドサットＴＭ画像上の任意の位置の座標（ピクセル番号とライン番号）であり、θ０は地図の座標系に対する画像の座標系の回転角であり、ΔはランドサットＴＭ画像の空間解像度である。
【００１１】
本発明では、数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、該シミュレーション画像及びランドサットＴＭ画像の画素値の相関係数を評価関数として、適切な平行移動パラメータを求める（評価関数の最適値が得られるパラメータを求める）。ここで太陽の直達入射照度画像とは、太陽の位置情報及び既知の標高情報から得た太陽の直達入射照度の差を含む画像である。太陽の直達入射照度ｃｏｓβは、ｃｏｓβ＝ｃｏｓθ＊ｃｏｓ（ｅ）＋ｓｉｎ（ｅ）＊ｃｏｓ（φ−Ａ）の式に基いて演算することができる。但しθは天頂角であり、Ａは方位角であり、ｅが地表面の斜度であり、φが傾斜方位である。画素値とは、画像中の１画素の明るさに比例するデジタル・ナンバーであり、ランドサットＴＭ画像の画素値は、１画素の輝度に比例する値であり、直達入射照度画像の画素値は、１画素の照度に比例する値である。
【００１２】
また相関係数は、ランドサットＴＭ画像及びシミュレーション画像中のそれぞれの画素値の関係性を評価する指数であり、−１から＋１の間の数値となる。相関係数は、一般的な公式に基いて求めることができる。そして相関係数を評価関数とするためには、まずランドサットＴＭ画像とシミュレーション画像の大まかなズレ量を目視等により求めておき、そのズレ量を中心にして一方の画像を他方の画像に対してピクセル方向及びライン方向に所定量ずつ移動させて、特定の画素領域（例えば画像の中心の１００×１００の画素の領域）における前述の相関係数を求める。そして、この相関係数が最大になる（最適値になる）移動量を決定して、この移動量から二つの画像のズレ量を求め、このズレ量を補正できるように平行移動パラメータの値を決定する。
【００１３】
直接的には、ｐ０及びｌ０を評価関数を用いて最適化または適切化するパラメータとすればよい。しかし具体的には、上記二つの式を、下記のように簡略化する。
【００１４】
ｐ＝ａｘ−ｂｙ＋ｃ
ｌ＝ｂｘ−ａｙ＋ｄ
但しａ乃至ｄは変数である。そしてこの式のｃ及びｄを平行移動パラメータとする。そしてこれらのパラメータを、二つの画像の画素値の相関係数を評価関数として、この評価関数が最適値になるように定め、前述のｐ０及びｌ０を求める。
【００１５】
このようにして定めた平行移動パラメータを用いて上記幾何変換式に従い、対象領域のすべての座標について幾何変換を行えば、多数のＧＣＰを取得することなく、対象領域全体にわたって高い精度の幾何補正を行うことができる。
【００１６】
なお平行移動パラメータをそれぞれ変化させたときのピクセル方向の相関係数及びライン方向の相関係数の値が最大値となる値を、平行移動パラメータの最適値とすれば、１画素内の微調整が可能になる。
【００１７】
本発明の特徴をより抽象的に表現すれば、本発明の特徴は、ランドサットＴＭ画像全体に現れる画像上の特徴（例えば、前述の輝度）と地図上の地図情報に関連して得られる地図関連情報のうち画像上の特徴に対応する対応地図情報（例えば、前述の太陽の直達入射照度）との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて評価を行いながら、平行移動パラメータの最適化を行うことである。ここで「空間的整合性を統計的に偏りなく定量的に評価する評価関数」の一つが、前述のＴＭ画像上の輝度に基く画素値と太陽の直達入射照度画像（シミュレーション画像）上の太陽の直達入射照度に基く画素値との間の相関係数である。このような評価関数としては、他のものを用いることもできる。例えば、ＴＭ画像上の特徴が、この画像上に現れる陸地と海との相違であり、対応地図情報が地図上に示される陸地と海との相違であるとした場合には、評価関数として、海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの識別データの正否を、対応地図情報により判定して得た識別率を用いることができる。尚「層別」とは、統計処理において、事前情報（海岸線からの距離）により、母集団をいくつかの部分集団（層またはカテゴリー）に分けることを意味する。識別率を評価関数とする場合にも、この評価関数が最適値または適切値になるようにパラメータを定める。
【００１８】
本発明の特徴を、衛星画像に対する精密幾何補正方法として更に抽象的に（一般的に）表現すれば、次のように表現できる。
【００１９】
すなわち本発明は、地図上の座標と衛星画像の座標との関係を表すために複数のパラメータを含んで構成された幾何変換式を用いて衛星画像の座標と地図の座標とを整合させるように衛星画像を幾何補正する衛星画像の精密幾何補正方法を対象として、衛星画像全体に現れる画像上の特徴と地図上の地図情報に関連して得られる地図関連情報のうち画像上の特徴に対応する対応地図情報との空間的整合性を統計的に偏りなく定量的に評価する評価関数を用いて複数のパラメータの最適化を行うことを特徴とするものである。
【００２０】
さらにパラメータの決定を容易にするためには、２つの段階を経てパラメータを決定すればよい。まず最初に、ランドサットＴＭ画像（衛星画像）上に海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの判定の正否を地図上に示される陸地と海のデータに基いて判定して、その識別率を求める。そしてこの識別率を評価関数として複数のパラメータの適切値を求める。次にこの適切値を幾何変換式に入力して衛星画像を幾何補正して幾何補正された衛星画像を定める。その後、数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、このシミュレーション画像及び幾何補正された衛星画像のそれぞれの画素値の相関係数を評価関数として、複数のパラメータの最適値化を行う。このようにすると目視によりランドサットＴＭ画像（衛星画像）とシミュレーション画像のズレを確認する必要がなくなるので、パラメータの決定の自動化が容易になる。
【００２１】
【発明の実施の形態】
以下図面を参照して本発明の実施の形態を詳細に説明する。図１は、本発明の第１の実施の形態の方法のデータ処理の流れを説明するための図である。この実施の形態では、衛星画像としてランドサットＴＭ画像を扱う。ランドサットＴＭデータのシステム情報に含まれる標準処理データに付属するリーダファイルの地図投影アンシナリーレコードに関する情報には、撮影シーンのセンター（中心）の基準となる地図上の中心座標即ちＵＴＭ座標（ｘ０，ｙ０）、ＵＴＭ座標系に対する送信するＴＭ画像の回転角（オリエンテーション角θ０）、および空間解像度（Δ）が含まれている。空間解像度は、画像を撮影したセンサが分解できる最も小さい対象物の単位である。これらの情報を用いれば、以下の幾何変換式（１）及び（２）を用いて任意の地図上のＵＴＭ座標ｘ，ｙを衛星画像上の位置（ラインおよびピクセルに）に変換（いわゆる逆変換）することができる。これらの式（１）及び（２）は、システム情報に付属している。
【００２２】
ｐ−ｐ０＝［ｃｏｓθ０＊（ｘ−ｘ０）−ｓｉｎθ０＊（ｙ−ｙ０）］／Δ  …（１）
ｌ−ｌ０＝［−ｓｉｎθ０＊（ｘ−ｘ０）−ｃｏｓθ０＊（ｙ−ｙ０）］／Δ   …（２）
上記式において、ｘ０及びｙ０は地図の座標系の中心座標であり、ｘ及びｙは地図の座標系の任意の位置の座標であり、ｐ０及びｌ０はランドサットＴＭ画像の座標系の中心座標（ピクセル番号とライン番号）であり、ｐ及びｌはランドサットＴＭ画像上の任意の位置の座標（ピクセル番号とライン番号）であり、θ０はオリエンテーション角であり、Δは空間解像度である。
【００２３】
上記式をアフィン変換式で簡略化すれば次のようになる。
【００２４】
ｐ＝ａｘ−ｂｙ＋ｃ
ｌ＝ｂｘ−ａｙ＋ｄ
但しａ乃至ｄは変数である。この場合ｃ及びｄが平行移動パラメータとなる。システム幾何変換では上式の係数はすべてが決定される。本実施の形態の方法では、画像座標上での平行移動を意味するｃとｄの２つのパラメータを変数と考える。また上記式（１）及び（２）におけるｐ０及びｌ０をパラメータ（変数）と考えても同じである。
【００２５】
この実施の形態では、ランドサットＴＭ画像を所定の幾何変換を経て幾何補正したもの（変換された衛星画像）とシミュレーションにより作成した直達入射照度画像との相関係数を幾何補正の評価関数として用いてパラメータの最適値を決定する。相関係数は位置決めが正しい程大きくなるため、精密幾何補正に用いる幾何変換式のパラメータを最適化することが可能となる。最適化は相関係数が大きくなるようにパラメータを探索することによって実現される。
【００２６】
なお、この実施の形態では、上記式を用いた幾何変換に正射投影を含めるため、上式で得られるピクセル値に起伏に基づく位置ズレを加えた衛星画像上の見かけの位置のデータを利用している。最適解を求めるために利用する幾何変換には正射投影を組み込むものとする。すなわち、上記（１）及び（２）式に地図座標（ｘ，ｙ）を代入して得られるピクセル値ｐに起伏に基づく位置ズレΔｐを加える。ただし、起伏に基づく位置ズレΔｐを厳密に求めるには、地球の曲率を考慮した計算が必要となり、計算機への負担が大きくなる。なお厳密な計算による位置ズレが地球の曲率を考慮しない場合の位置ズレよりも１０％程大きくなることを利用することにより計算時間の短縮をはかることができる。また衛星直下点の起伏に基く位置ズレは、原画像における欠損データの軌跡から推定すればよい。なおこの点に関しては、「写真測量とリモートセンシング」のＶｏｌ．３７，Ｎｏ．４の１２頁〜２２頁に掲載された「ランドサットＴＭ画像の正射投影とその評価」と題する論文に掲載されている。
【００２７】
シミュレーション画像を得るために必要な太陽の直達入射照度ｃｏｓβは、太陽位置（天頂角θと方位角Ａ）と地表面の斜度ｅと傾斜方位φの関数として以下の式から与えられる。
【００２８】
ｃｏｓβ＝ｃｏｓθｃｏｓ（ｅ）＋ｓｉｎθｓｉｎ（ｅ）ｃｏｓ（φ−Ａ）
ここで太陽位置（天頂角θと方位角Ａ）のデータについては、標準処理データに付属する太陽位置を利用する。斜度ｅと傾斜方位φを求めるには数値標高モデルが必要となる。数値標高モデルは、国土地理院から発行されている数値地図５０ｍ（標高）を投影変換（ＵＴＭ座標系）および再配列（３０ｍ）して用いる。国土地理院から発行されている数値地図５０ｍメッシュ（標高）は、格子点（メッシュ）が等間隔の緯度と経度により構成されており、サイズが約５０ｍである。したがって、これをランドサットＴＭの画像処理に利用するには、投影変換（ＵＴＭ座標系）および再配列（空間解像度３０ｍ）が必要になるのである。なおこの点の詳細に関しては、日本リモートセンシング学会誌、Ｖｏｌ．２１，Ｎｏ．２の１５０〜１５７頁に掲載された「数値標高モデルの投影変換に用いる内挿方の評価法」と題する論文に掲載されている。本実施の形態では、この投影変換（ＵＴＭ座標系）および再配列（空間解像度３０m）したものから、斜度ｅと傾斜方位φを求めている。
【００２９】
本実施の形態では、図２に示す岩手県内の５箇所（黒丸の中に１から５の数字を記載した地域）を対象地域とした。各地点毎に切り出す衛星画像の大きさは縦４３０、横４３０、空間解像度は３０ｍとして数値標高モデルを作成した。
【００３０】
衛星画像は１９９５年６月１２日撮影のランドサットＴＭデータ（パス１０７，ロウ３２）を用いた。撮影時の太陽高度は５８度、方位は北から時計周りの方向に１１２度であった。この情報を用いて太陽の直達入射照度のシミュレーション画像を作成した。地点１のシミュレーション画像を図３に示す。
【００３１】
最初にシステム幾何変換後の衛星画像とシミュレーション画像とを目視で比較することにより、おおまかな位置合わせを行った。その結果ピクセル方向に４０画素、ライン方向に２７画素の位置ズレが確認された。そこで目視により確認した位置ズレ（ピクセル方向に４０画素、ライン方向に２７画素）を中心にして、変換後のランドサットＴＭ画像をシミュレーション画像に対して画素単位で（１画素ずつ）ピクセル方向及びライン方向に移動させて、ランドサットＴＭ画像の中心に位置する例えば１００×１００個の画素の画素値とそれらに対応するシュミレーション画像の画素の画素値との相関係数を求めた。その結果を、下記表１に示す。
【００３２】
【表１】
0004005336.tif 000002
上記表１の横軸の数字はピクセル方向に４０画素を中心にして±２画素ずらすことを意味し、縦軸の数字はライン方向に２７画素を中心にして±２画素ずらすことを意味する。そして表中の数字は、画素単位で移動したときの画素値の相関係数である。表１の結果が、評価関数を示している。表１からは、ピクセル方向に４０画素、ライン方向に２７画素ずれた位置における相関係数０．６９５が最も大きくなることが分かる。すなわち評価関数の最適値は、ピクセル方向に４０画素、ライン方向に２７画素ずらすことである。この段階で、このようなズレ量を得られるように前述の平行移動パラメータｃ及びｄを決定することができる。
【００３３】
図４には地点１の変換後の衛星画像（バンド５）を示した。図３に示したシミュレーション画像と比べて、両者の陰影が一致していることが確認できる。
【００３４】
次に、１画素以内の微調整を行うために行う、前述の位置ズレを原点とした幾何変換式からの平行移動パラメータの最適値の決定方法を説明する。平行移動パラメータｃ，ｄを変化させて作成した衛星画像（ピクセル方向に４０画素、ライン方向に２７画素ずらしたランドサットＴＭ画像）の中心領域の１００×１００の画素の画素値とこれに対応するシミュレーション画像の画素の画素値の相関係数を計算した。その結果の例を図５及び図６に示す。ピクセル方向への移動量を決める平行移動パラメータｃを変化させた場合の相関係数の変化を図５に示しており、ライン方向への移動量を決める平行移動パラメータｄを変化させた場合の相関係数の変化を図６に示してある。図５及び図６を比較すると、パラメータｃを変化させた場合のほうが、パラメータｄを変化させた場合よりも変化がより明らかであるが、両者ともきれいな対称性がある。２つのパラメータを２次元的に探索することにより最適なパラメータを１画素よりも一桁小さいオーダーで決定することができる。これは相関係数の計算に非常に多くの点を利用できるため、結果の精度が向上したためであると思われる。すなわちパラメータｃ及びｄを変化させて相関係数が最大値になるときのそれぞれのパラメータｃ及びｄの値（この例では０．２及び０．２）をそれぞれの最適値として定めることができる。なお、相関係数を計算する画像範囲を地上被覆が均一で起伏の比較的大きな範囲とすれば０．８をこえる相関係数が得られる。しかしながら、そのようにしても位置ズレの情報としては大きな変化は見られないため、単純に標本の数を大きくするという意味で相関係数をとる範囲は全域にとるのが好ましい。
【００３５】
このようにして最適な平行移動量を図１の各地点１乃至５に対して求めた結果を表２に示す。
【００３６】
【表２】
0004005336.tif 000003
表２において、横軸のｃ及びｄはそれぞれパラメータを示し、ｒは相関係数の値を示している。選択した５点はかなり広い範囲を占めるが、最適な平行移動パラメータｃ及びｄの違いはピクセル方向、ライン方向ともで１．１画素程度であった。この程度の範囲ではｃとｄを固定しても１画素以内の、ライン方向で約１画素に満たない。これはシステム情報に基づく幾何補正の精度がかなり良いことを示していると考えられる。また、西側の地点（図１の地点２と５）でピクセル方向のマイナス側へのズレが大きくなり、ライン方向ではプラス側へのズレが大きくなるなど系統的な特徴が見られる。このことは、より大きな範囲を問題にする場合には回転を含めた最適化が必要なことを示唆している。
【００３７】
上記の説明では、パラメータｃ及びｄの最適値を求めた。しかし上記式（１）及び（２）中のｐ０及びｌ０即ちシーンのセンターをパラメータとすることができる。この場合には、幾何変換式のシーンのセンター（ｐ０、ｌ０）を変化(δｐ、δｌ)させることにより１画素より小さな単位で位置ズレを調整する。また、正射投影を行うことにより、広い領域で精度の高い出力画像を作成することが期待できる。更に最適解の計算は、ＩＤＬ（Interactive Data Language)（Ver5.3)に組み込まれているPowellの方法を用いて行うことができる。
【００３８】
図１の例では、パラメータの最適値化を行った後に式（１）及び（２）の幾何変換式にパラメータを代入して、ランドサットＴＭ画像の精密幾何変換を行っている。しかしながら図７に示すように、ランドサットＴＭデータを幾何変換する際に、式（１）及び（２）を用い、最初に適当な平行移動パラメータの値をこれらの式中に挿入して、パラメータの最適値化を実行し、その結果で平行移動パラメータを変更するようにしてもよい。
【００３９】
本実施の形態では、ランドサットＴＭの標準処理データを幾何変換する簡便な方法（修正システム幾何変換法）が提案されている。この修正システム幾何変換法ではシステム情報に含まれる地図投影に関する情報を基本に、平行移動と地上の起伏による位置ズレ（正射投影）を加味した変換を考えている。そして、太陽の直達入射照度のシミュレーション画像と衛星画像の画素値の相関係数を評価関数として平行移動のパラメータを調整する。０．１画素単位の微調整も可能である。従来のＧＣＰを利用した精密幾何補正では、ＧＣＰの選択と座標変換式の統計的な推定という２段階の操作が必要となるのに対して、本発明の方法では従来問題とされてきたＧＣＰの選択が省略できる利点がある。
【００４０】
上記実施の形態では、２つの平行移動パラメータの最適値化を行ったが、最適値化できるパラメータは平行移動パラメータに限定されるものではなく、他のパラメータについても同様の方法で最適値化を実行できるのは勿論である。
【００４１】
図８は、本発明の第２の実施の形態のデータ処理の流れを示す図である。上記の第１の実施の形態では、目視で衛星画像とシミュレーション画像とのズレ量を見てから、パラメータの最適値化を行っている。この実施の形態では、目視によるズレ量の検出に変えて、自動的に大域的なズレ量を検出する。そこでこの実施の形態では、大域的な方法（海と陸の違いを利用した方法）で、大域的なズレ量を検出した後、前述の第１の実施の形態で採用している局所的な方法（太陽の直達入射照度の推定値を利用した方法）でパラメータを決定する。
【００４２】
具体的に説明すると、この実施の形態では、最初にランドサットＴＭ画像上の海岸線からの距離で層別したカテゴリー毎に無作為に抽出した複数の評価点が陸地であるか海であるかの判定の正否を地図上に示される陸地と海のデータに基いて判定して、その識別率を求める。そしてこの識別率を評価関数として平行移動パラメータの適切値を求める。その後適切値を幾何変換式に入力してランドサットＴＭ画像を幾何補正して幾何補正ランドサットＴＭ画像を定める。次に数値標高モデルを用いて作成した太陽の直達人射照度画像をシミュレーション画像とし、シミュレーション画像及び幾何補正ランドサットＴＭ画像のそれぞれの画素値の相関係数を評価関数として、平行移動パラメータの最適値化を行う。
【００４３】
以下に大域的な方法（海と陸の違いを利用した方法）について説明する。衛星画像で識別できるもっとも大きな特徴は、陸域と海域の違いであり、これらは特に近赤外バンドの画像では明瞭に識別できる。地図座標上の各点が海か陸かの判定は国土地理院から発行される数値情報によって知ることができる。この特徴を利用するため、陸域と海域から複数の評価点を選択する。次に設定した幾何補正式を利用して評価点の衛星画像上の座標を求める。この座標上の輝度値から判定できる海と陸が地図情報と一致する評価点の割り合いを適合度とする。この方法では、評価点の設定が重要である。海岸に近いところだけでなく遠いところも含めることによりパラメータの初期値によらない最適化が可能となる。また評価点を調整することにより、計算時間の短縮を計ることができるため、シーン全体などの広い地域を対象にする場合や多数のパラメータを推定する場合に有効である。
【００４４】
この実施の形態では、宇宙開発事業団から提供されたランドサット５号からの衛星画像（ＴＭデータ）「１９９５年６月１２日撮影（パス１０７／ロウ３２）」を用いた。ここで使用するデータはいわゆる標準処理データ（レベル１ｂ）といわれるもので、（１）検知素子間の感度差の調整、および（２）衛星の軌道・姿勢情報を用いたＵＴＭ座標系への投影と再配列（システム幾何補正）が行われたものである。しかし、システム幾何補正では（１）地表面の起伏の影響が考慮されていない、（２）軌道情報が正確ではない、（３）画像の配列がＵＴＭ座標系のｘｙ軸と対応していない（オリエンテーション角）などの問題があるため、地図のデータと正確には重ね合わせることができない。
【００４５】
図９に可視バンドの画像ファイルの一例の画像を示している。ファイルの最初にはヘッダレコード（１行）、各レコードの最初の３２バイトはレコードヘッダ、最後の６８バイトはレコードトレイラであるが、本解析に必要な情報は含まれていない。他の部分が画像データであるが、左右の黒い部分（輝度値０）はシステム補正の際に生じた撮影されていない部分である。精密幾何変換（正射投影）では画像を各ライン毎のこの境界の軌跡の中間点を衛星直下点とする。
【００４６】
標準処理データには、画像ファイルの他に情報ファイル（リードファイル）が付属している。このファイルには幾何変換に必要な以下の情報が記載されている。
【００４７】
太陽高度  θ＝ 62 度
太陽方位  φ＝ 119 度
オリエンテーション角θ０＝ 0.1830542 radian
シーン中心（ＵＴＭ座標系）
ｘ０＝ 534146.1182 m
ｙ０＝ 4463628.9062 m
空間分解能  ２８．５ｍ
また、この画像ファイル上のシーン中心は
ｐ０＝ 3492.0 pixel
ｌ０＝ 2983.5 line
である。
【００４８】
これらの情報からＵＴＭ座標（ｘ、ｙ）が与えられた場合の画像上の位置（ｐ、ｌ）は上記第１の実施の形態で用いた（１）式及び（２）式を用いて次の式で求められる。
【００４９】
ｐ＝ｐ０＋［ｃｏｓθ０＊（ｘ−ｘ０）−ｓｉｎθ０＊（ｙ−ｙ０）］／28.5
ｌ＝ｌ０＋［−ｓｉｎθ０＊（ｘ−ｘ０）−ｃｏｓθ０＊（ｙ−ｙ０）］／28.5
図１０に地図上に衛星データの範囲を示した。図１０において、領域Ａが図９の衛星画像の範囲であり、領域Ｂが大地的な方法を実施する範囲であり、領域Ｃが局所的な方法を実施する範囲である。
【００５０】
図８に示したデータ処理の流れは、（Ａ）評価用のデータの作成と（Ｂ）評価関数の最適化による幾何変換式の推定に分けられる。後者には海と陸の判別を用いる大域的な方法と直達入射照度を用いる局所的な方法がある。それぞれ独立に用いることもできるが、ランドサットＴＭの場合には、大域的な方法で回転を含む幾何変換式を作成した後で、その平行移動パラメータを修正する形で局所的な方法を用いるのが合理的である。
【００５１】
図８のデータ処理の流れでは、まず最初に、評価用のデータを作成する。評価用のデータは２種類（海と陸の判別データと太陽の直達入射照度）あるが、どちらも国土地理院発行の数値地図５０mメッシュ（標高）から作成する。
【００５２】
図１１に示した対象とする衛星画像（パス１０７／ロウ３２用）の範囲と海岸線の多くが重なる矩形の領域Ｂの数値標高モデルを切り出したものが図１１に示すものである。この画像を用いて海岸線からの距離画像（海と陸が接するものをレベル１）を作成した。海岸線からの距離が約７５０ｍの領域を距離で１５のレベルに層別し、各レベル（カテゴリー）毎に２００画素（海と陸が半々）を無作為に抽出し、合計３０００個の地点を評価点とし表３に示すようなデータを作成した。この一部を示したものが図１２である。図１２において、＋印が評価点を示している。さらに幾何変換式の決定で用いるデータはこの経緯度で示される位置データをＵＴＭ座標系に変換する。
【００５３】
【表３】
0004005336.tif 000004
図１３は、図１０に示した地域Ｃの太陽の直達入射照度画像（シミュレーション画像）を示している。この画像は、数値標高モデル（ＵＴＭ座標系／空間分解能３０ m）を数値地図５０mメッシュ（標高）を投影変換／再配列（共１次内挿法）して作成したものである。画像の大きさは４３０×４３０である。太陽の直達入射照度ｃｏｓβの求め方は第1の実施の形態の場合と同じである。
【００５４】
次に評価関数の最適化による幾何変換式の推定について説明する。
【００５５】
図９に示したバンド４の衛星画像に対して、大域的な方法を適用する。システム情報に基づいて評価点の衛星画像上の位置を計算し、各評価点の輝度値（画素値）が閾値γ（２０を設定）以下の場合を海域とし、輝度値が閾値より大きい場合を陸域と判断した。そしてこの判断結果（識別データ）の正否を、前述の数値地図５０ｍメッシュ（標高）から得た地図関連情報に照らして判断し、識別率を求めた。その結果、評価点から衛星画像の範囲外にある１２２地点を除いた２８７８地点の中で７１８地点が間違って判定された。
【００５６】
この例では、間違った数を識別率とし、これを評価関数として評価値が大きくなるように（識別率が小さくなるように）シーンのセンターのズレ（Δｐ、Δｌ）および閾値（γ）に関して最適化を行ったところ、間違いは５１個に減少した。ここで最適値化は、Δｐ、Δｌおよび閾値γを周知の最適化手法に従って可変し、その都度、輝度値（画素値）と閾値との対比して行った。このときのパラメータの値は以下の通りである。
【００５７】
ケースＡ：  Δｐ＝40.90、Δｌ＝27.18、γ＝28.39
ケースＡは平行移動のみであるが、回転を含めるために、さらに、オリエンテーション角のズレΔθを含めて最適化を行ったところ、間違いは４０個となった。図１４に間違った場所と範囲外の地点を地図上に示した。北部に見えるのは衛星画像の範囲外となった評価点である。このときのパラメータの値は以下の通りである。
【００５８】
ケースＢ：  Δｐ＝40.25、Δｌ＝27.18、Δθ＝0.000217、γ＝29.08
空間分解能２８．５ｍも変化させて最適化を行った場合には、間違いの現象はみられなかった。なお、間違った地点を個別的にチェックはしていないが、陸奥小川原湖周辺の湖沼の影響が多く見られるほかは、系統的な傾向は見られない。小さな島や河口、あるいは港の回収などの影響によると思われる。
【００５９】
前述の大域的な方法で得られた幾何変換式を用いて、画像の再配列を行った。すなわち幾何補正したランドサットＴＭ画像とシミュレーション画像との間の大まかなズレ量を決定した。この画像の再配列が、第１の実施の形態における目視によるシーンのセンターのズレ量の決定に相当する。その後第１の実施の形態と同様にして、画素単位で幾何補正したランドサットＴＭ画像をピクセル方向及びライン方向に移動して、変換された衛星画像とシミュレーション画像の所定領域におけるそれぞれの画素値の相関係数の変化を求めた。表４及び表５は、それぞれ前述のケースＡ及びＢのパラメータを用いた場合に得られた相関係数の変化を示している。なおこの例では、地形の起伏による誤差を避けるために正射投影を考慮した変換を行った。また内挿には共１次内挿法を用いた。解析には変換後の衛星画像をピクセル方向とライン方向にずらして求めた相関係数を用いた。正しい変換であればズレがないところで最大値をとる。
【００６０】
【表４】
0004005336.tif 000005
【表５】
0004005336.tif 000006
表４および表５には、ともに横軸にピクセル方向への移動量（画素数）を示して縦軸にライン方向への移動量（画素数）を示している。横軸及び縦軸の数字は、ケースＡ及びケースＢで求めたシーンのセンターのズレΔｐ及びΔｌの値を基準にして（０にして）画素単位で画像をピクセル方向及びライン方向にずらすことを意味している。
【００６１】
ケースＡ（平行移動のみ）の場合も、ケースＢ（回転も含める）の場合も、共にライン方向にはズレがなく、ピクセル方向に最大２画素程度のズレが見られる。両者を比較すると、後者の方の位置ズレが若干小さい。
【００６２】
次に相関係数を評価関数として、各ケースについてパラメータの最適化を行った。その結果を以下に示す。
【００６３】
ケースＡ：  Δｐ＝38.79、Δｌ＝27.45、相関係数0.692
ケースＢ：  Δｐ＝39.19、Δｌ＝27.21、Δθ＝−0.000213、相関係数0.692
最適化を行った場合には、回転を含める効果はほとんど現れない。これは対象画像の範囲（約１０ｋｍ）では回転の影響は数ｍ（１００００×Δθ）程度しかあらわれないためである。したがって、回転のズレを求めるには大域的な方法が不可欠である。また、局所的な方法で最適化を行う場合、初期値を適切に設定しないと解がなかなか見つからないことが起きる。大域的な方法で得られたパラメータを初期値として利用することにより、この問題を解決することができる。
【００６４】
【発明の効果】
本発明によれば、相関関数を評価関数として、適切な平行移動パラメータを求めることができるので、この平行移動パラメータを用いて対象領域のすべての座標について幾何変換を行えば、多数のＧＣＰを取得することなく、対象領域全体にわたって高い精度の幾何補正を行うことができる利点がある。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態の方法のデータ処理の流れを説明するための図である。
【図２】対象地域を示す図である。
【図３】図２に示した地点１のシミュレーション画像を示す図である。
【図４】図２に示した地点１の変換後の衛星画像（バンド５）を示す図である。
【図５】パラメータｃに対する相関係数の変化を示す図である。
【図６】パラメータｄに対する相関係数の変化を示す図である。
【図７】第１の実施の形態の別のデータ処理の流れを示す図である。
【図８】本発明の第２の実施の形態のデータ処理の流れを示す図である。
【図９】可視バンドの画像ファイルの画像を示す図である。
【図１０】対象地域を示す図である。
【図１１】図１０の領域Ｂの数値標高モデルを示す図である。
【図１２】無作為に抽出した評価点の例を示す図である。
【図１３】図１０に示した地域Ｃの直達入射照度画像（シミュレーション画像）を示す図である。
【図１４】陸と海の判定を誤った評価点を示す図である。
【符号の説明】
Ａ乃至Ｃ  領域
【図１】
0004005336.tif 000007
【図２】
0004005336.tif 000008
【図３】
0004005336.tif 000009
【図４】
0004005336.tif 000010
【図５】
0004005336.tif 000011
【図６】
0004005336.tif 000012
【図７】
0004005336.tif 000013
【図８】
0004005336.tif 000014
【図９】
0004005336.tif 000015
【図１０】
0004005336.tif 000016
【図１１】
0004005336.tif 000017
【図１２】
0004005336.tif 000018
【図１３】
0004005336.tif 000019
【図１４】
0004005336.tif 000020