JP4005336B2 - Precision geometric correction method for Landsat TM image and precise geometric correction method for satellite image - Google Patents

Description

上記表１の横軸の数字はピクセル方向に４０画素を中心にして±２画素ずらすことを意味し、縦軸の数字はライン方向に２７画素を中心にして±２画素ずらすことを意味する。そして表中の数字は、画素単位で移動したときの画素値の相関係数である。表１の結果が、評価関数を示している。表１からは、ピクセル方向に４０画素、ライン方向に２７画素ずれた位置における相関係数０．６９５が最も大きくなることが分かる。すなわち評価関数の最適値は、ピクセル方向に４０画素、ライン方向に２７画素ずらすことである。この段階で、このようなズレ量を得られるように前述の平行移動パラメータｃ及びｄを決定することができる。
【００３３】
図４には地点１の変換後の衛星画像（バンド５）を示した。図３に示したシミュレーション画像と比べて、両者の陰影が一致していることが確認できる。
【００３４】
次に、１画素以内の微調整を行うために行う、前述の位置ズレを原点とした幾何変換式からの平行移動パラメータの最適値の決定方法を説明する。平行移動パラメータｃ，ｄを変化させて作成した衛星画像（ピクセル方向に４０画素、ライン方向に２７画素ずらしたランドサットＴＭ画像）の中心領域の１００×１００の画素の画素値とこれに対応するシミュレーション画像の画素の画素値の相関係数を計算した。その結果の例を図５及び図６に示す。ピクセル方向への移動量を決める平行移動パラメータｃを変化させた場合の相関係数の変化を図５に示しており、ライン方向への移動量を決める平行移動パラメータｄを変化させた場合の相関係数の変化を図６に示してある。図５及び図６を比較すると、パラメータｃを変化させた場合のほうが、パラメータｄを変化させた場合よりも変化がより明らかであるが、両者ともきれいな対称性がある。２つのパラメータを２次元的に探索することにより最適なパラメータを１画素よりも一桁小さいオーダーで決定することができる。これは相関係数の計算に非常に多くの点を利用できるため、結果の精度が向上したためであると思われる。すなわちパラメータｃ及びｄを変化させて相関係数が最大値になるときのそれぞれのパラメータｃ及びｄの値（この例では０．２及び０．２）をそれぞれの最適値として定めることができる。なお、相関係数を計算する画像範囲を地上被覆が均一で起伏の比較的大きな範囲とすれば０．８をこえる相関係数が得られる。しかしながら、そのようにしても位置ズレの情報としては大きな変化は見られないため、単純に標本の数を大きくするという意味で相関係数をとる範囲は全域にとるのが好ましい。
【００３５】
このようにして最適な平行移動量を図１の各地点１乃至５に対して求めた結果を表２に示す。
【００３６】
【表２】

表２において、横軸のｃ及びｄはそれぞれパラメータを示し、ｒは相関係数の値を示している。選択した５点はかなり広い範囲を占めるが、最適な平行移動パラメータｃ及びｄの違いはピクセル方向、ライン方向ともで１．１画素程度であった。この程度の範囲ではｃとｄを固定しても１画素以内の、ライン方向で約１画素に満たない。これはシステム情報に基づく幾何補正の精度がかなり良いことを示していると考えられる。また、西側の地点（図１の地点２と５）でピクセル方向のマイナス側へのズレが大きくなり、ライン方向ではプラス側へのズレが大きくなるなど系統的な特徴が見られる。このことは、より大きな範囲を問題にする場合には回転を含めた最適化が必要なことを示唆している。
【００３７】
上記の説明では、パラメータｃ及びｄの最適値を求めた。しかし上記式（１）及び（２）中のｐ_０及びｌ_０即ちシーンのセンターをパラメータとすることができる。この場合には、幾何変換式のシーンのセンター（ｐ_０、ｌ_０）を変化(δｐ、δｌ)させることにより１画素より小さな単位で位置ズレを調整する。また、正射投影を行うことにより、広い領域で精度の高い出力画像を作成することが期待できる。更に最適解の計算は、ＩＤＬ（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メッシュ（標高）から作成する。
【００５２】
図１１に示した対象とする衛星画像（パス１０７／ロウ３２用）の範囲と海岸線の多くが重なる矩形の領域Ｂの数値標高モデルを切り出したものが図１１に示すものである。この画像を用いて海岸線からの距離画像（海と陸が接するものをレベル１）を作成した。海岸線からの距離が約７５０ｍの領域を距離で１５のレベルに層別し、各レベル（カテゴリー）毎に２００画素（海と陸が半々）を無作為に抽出し、合計３０００個の地点を評価点とし表３に示すようなデータを作成した。この一部を示したものが図１２である。図１２において、＋印が評価点を示している。さらに幾何変換式の決定で用いるデータはこの経緯度で示される位置データをＵＴＭ座標系に変換する。
【００５３】
【表３】

図１３は、図１０に示した地域Ｃの太陽の直達入射照度画像（シミュレーション画像）を示している。この画像は、数値標高モデル（ＵＴＭ座標系／空間分解能３０ m）を数値地図５０mメッシュ（標高）を投影変換／再配列（共１次内挿法）して作成したものである。画像の大きさは４３０×４３０である。太陽の直達入射照度ｃｏｓβの求め方は第1の実施の形態の場合と同じである。
【００５４】
次に評価関数の最適化による幾何変換式の推定について説明する。
【００５５】
図９に示したバンド４の衛星画像に対して、大域的な方法を適用する。システム情報に基づいて評価点の衛星画像上の位置を計算し、各評価点の輝度値（画素値）が閾値γ（２０を設定）以下の場合を海域とし、輝度値が閾値より大きい場合を陸域と判断した。そしてこの判断結果（識別データ）の正否を、前述の数値地図５０ｍメッシュ（標高）から得た地図関連情報に照らして判断し、識別率を求めた。その結果、評価点から衛星画像の範囲外にある１２２地点を除いた２８７８地点の中で７１８地点が間違って判定された。
【００５６】
この例では、間違った数を識別率とし、これを評価関数として評価値が大きくなるように（識別率が小さくなるように）シーンのセンターのズレ（Δｐ、Δｌ）および閾値（γ）に関して最適化を行ったところ、間違いは５１個に減少した。ここで最適値化は、Δｐ、Δｌおよび閾値γを周知の最適化手法に従って可変し、その都度、輝度値（画素値）と閾値との対比して行った。このときのパラメータの値は以下の通りである。
【００５７】
ケースＡ： Δｐ＝40.90、Δｌ＝27.18、γ＝28.39
ケースＡは平行移動のみであるが、回転を含めるために、さらに、オリエンテーション角のズレΔθを含めて最適化を行ったところ、間違いは４０個となった。図１４に間違った場所と範囲外の地点を地図上に示した。北部に見えるのは衛星画像の範囲外となった評価点である。このときのパラメータの値は以下の通りである。
【００５８】
ケースＢ： Δｐ＝40.25、Δｌ＝27.18、Δθ＝0.000217、γ＝29.08
空間分解能２８．５ｍも変化させて最適化を行った場合には、間違いの現象はみられなかった。なお、間違った地点を個別的にチェックはしていないが、陸奥小川原湖周辺の湖沼の影響が多く見られるほかは、系統的な傾向は見られない。小さな島や河口、あるいは港の回収などの影響によると思われる。
【００５９】
前述の大域的な方法で得られた幾何変換式を用いて、画像の再配列を行った。すなわち幾何補正したランドサットＴＭ画像とシミュレーション画像との間の大まかなズレ量を決定した。この画像の再配列が、第１の実施の形態における目視によるシーンのセンターのズレ量の決定に相当する。その後第１の実施の形態と同様にして、画素単位で幾何補正したランドサットＴＭ画像をピクセル方向及びライン方向に移動して、変換された衛星画像とシミュレーション画像の所定領域におけるそれぞれの画素値の相関係数の変化を求めた。表４及び表５は、それぞれ前述のケースＡ及びＢのパラメータを用いた場合に得られた相関係数の変化を示している。なおこの例では、地形の起伏による誤差を避けるために正射投影を考慮した変換を行った。また内挿には共１次内挿法を用いた。解析には変換後の衛星画像をピクセル方向とライン方向にずらして求めた相関係数を用いた。正しい変換であればズレがないところで最大値をとる。
【００６０】
【表４】

【表５】

表４および表５には、ともに横軸にピクセル方向への移動量（画素数）を示して縦軸にライン方向への移動量（画素数）を示している。横軸及び縦軸の数字は、ケースＡ及びケースＢで求めたシーンのセンターのズレΔｐ及びΔｌの値を基準にして（０にして）画素単位で画像をピクセル方向及びライン方向にずらすことを意味している。
【００６１】
ケースＡ（平行移動のみ）の場合も、ケースＢ（回転も含める）の場合も、共にライン方向にはズレがなく、ピクセル方向に最大２画素程度のズレが見られる。両者を比較すると、後者の方の位置ズレが若干小さい。
【００６２】
次に相関係数を評価関数として、各ケースについてパラメータの最適化を行った。その結果を以下に示す。
【００６３】
ケースＡ： Δｐ＝38.79、Δｌ＝27.45、相関係数0.692
ケースＢ： Δｐ＝39.19、Δｌ＝27.21、Δθ＝−0.000213、相関係数0.692
最適化を行った場合には、回転を含める効果はほとんど現れない。これは対象画像の範囲（約１０ｋｍ）では回転の影響は数ｍ（１００００×Δθ）程度しかあらわれないためである。したがって、回転のズレを求めるには大域的な方法が不可欠である。また、局所的な方法で最適化を行う場合、初期値を適切に設定しないと解がなかなか見つからないことが起きる。大域的な方法で得られたパラメータを初期値として利用することにより、この問題を解決することができる。
【００６４】
【発明の効果】
本発明によれば、相関関数を評価関数として、適切な平行移動パラメータを求めることができるので、この平行移動パラメータを用いて対象領域のすべての座標について幾何変換を行えば、多数のＧＣＰを取得することなく、対象領域全体にわたって高い精度の幾何補正を行うことができる利点がある。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態の方法のデータ処理の流れを説明するための図である。
【図２】対象地域を示す図である。
【図３】図２に示した地点１のシミュレーション画像を示す図である。
【図４】図２に示した地点１の変換後の衛星画像（バンド５）を示す図である。
【図５】パラメータｃに対する相関係数の変化を示す図である。
【図６】パラメータｄに対する相関係数の変化を示す図である。
【図７】第１の実施の形態の別のデータ処理の流れを示す図である。
【図８】本発明の第２の実施の形態のデータ処理の流れを示す図である。
【図９】可視バンドの画像ファイルの画像を示す図である。
【図１０】対象地域を示す図である。
【図１１】図１０の領域Ｂの数値標高モデルを示す図である。
【図１２】無作為に抽出した評価点の例を示す図である。
【図１３】図１０に示した地域Ｃの直達入射照度画像（シミュレーション画像）を示す図である。
【図１４】陸と海の判定を誤った評価点を示す図である。
【符号の説明】
Ａ乃至Ｃ 領域[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a precise geometric correction method for Landsat TM images and a precise geometric correction method for satellite images.
[0002]
[Prior art]
The standard processing data of Landsat TM data (Thematic Mapper Data) includes, for example, a geometric transformation formula including a plurality of parameters as information (system information) for projecting a satellite image on a map (UTM coordinate system). . However, in the geometric transformation based on the system information (system geometric transformation), since a considerable error occurs, it is generally rarely used. The Landsat TM data available from the Space Development Agency is an image obtained by performing a predetermined geometric correction on an actual satellite image. However, the geometric correction accuracy of the converted image is not necessarily high. Therefore, conventionally, when geometric accuracy is required, a precise geometric correction method (including orthographic projection) based on affine transformation using a plurality of GCPs (ground control points) is generally employed. is there. In this precise geometric correction method using GCP (ground control points), an affine transformation equation is statistically calculated by using the least square method from coordinates on a map corresponding to satellite images of a plurality of GCPs identified by visual observation. To estimate. Since an error cannot be avoided in determining the GCP for each pixel, it has been considered that it is necessary to acquire as many GCPs as possible in order to statistically improve the precision of the precise geometric correction. However, GCP acquisition is inefficient and there is a problem that high accuracy cannot be expected unless it is performed by a skilled person.
[0003]
Therefore, the inventors have used a solar direct irradiance image that can be created using a digital elevation model as a method for automatically determining a large number of GCPs as a template (simulation image), and a satellite image after system geometric transformation. A method of local matching (simulation image method) was proposed ("Photogrammetry and remote sensing" Vol. 37, No. 6, 1998 entitled "Identification of ground control points in satellite images using a digital elevation model") Announced).
[0004]
[Problems to be solved by the invention]
As a result of applying the simulation image method previously proposed by the inventor to a number of Landsat TM images, if this method is adopted, if the target area is not so wide, the size of the positional deviation of the satellite image by system geometric transformation It turns out that the direction is almost independent of location. This is presumably because Landsat equipped with a star sensor is controlled with extremely high accuracy.
[0005]
However, in this conventional simulation image method, satellite images can be locally projected (matched) on the map (the coordinate system of the satellite image and the coordinate system of the map can be locally matched), but the satellite in the entire target area The image could not be projected onto the map.
[0006]
An object of the present invention is to provide a Landsat TM image accurate geometric correction method and a satellite image accurate geometric correction method capable of performing accurate geometric correction on the entire target area without acquiring a large number of GCPs. There is to do.
[0007]
[Means for Solving the Problems]
In the present invention, the parameters of the geometric transformation formula used in the system geometric transformation are optimized. At a specific level, the translation parameter is optimized. At this time, parameters are selected (optimized) using the correlation coefficient between the satellite image and the simulation image as an evaluation function. Since the parameters to be changed are only two parallel movement parameters, it is easy to search for an optimal solution and fine adjustment within one pixel is possible.
[0008]
Therefore, the present invention matches the coordinates of the Landsat TM image and the coordinates of the map by using a geometric transformation formula including a plurality of parameters in order to express the relationship between the coordinates on the map and the coordinates of the Landsat TM image. Therefore, an object of improvement is a precise geometric correction method for Landsat TM images in which the Landsat TM images are geometrically corrected. In the present invention, a geometric transformation formula including a translation parameter included in system information provided along with the Landsat TM image is used as the geometric transformation formula.
[0009]
For example, the geometric transformation formula included in the system information is the following formula.
[0010]
p-p ₀ = [Cos θ ₀ * (Xx ₀ ) -Sinθ ₀ * (Y-y ₀ )] / Δ
l-1 ₀ = [-Sinθ ₀ * (Xx ₀ ) -Cosθ ₀ * (Y-y ₀ )] / Δ
In the above formula, x ₀ And y ₀ Is the center coordinate of the map coordinate system, x and y are the coordinates of any position in the map coordinate system, and p ₀ And l ₀ Is the center coordinates (pixel number and line number) of the coordinate system of the Landsat TM image, p and l are coordinates (pixel number and line number) of an arbitrary position on the Landsat TM image, and θ ₀ Is the rotation angle of the image coordinate system relative to the map coordinate system, and Δ is the spatial resolution of the Landsat TM image.
[0011]
In the present invention, an appropriate parallel movement parameter is obtained by using a solar direct irradiance image created using a digital elevation model as a simulation image, and using a correlation coefficient between pixel values of the simulation image and the Landsat TM image as an evaluation function. (A parameter for obtaining the optimum value of the evaluation function is obtained). Here, the direct incident illuminance image of the sun is an image including a difference in direct incident illuminance of the sun obtained from the position information of the sun and the known altitude information. The direct incident illuminance cosβ of the sun can be calculated based on the equation: cosβ = cosθ * cos (e) + sin (e) * cos (φ−A). Where θ is the zenith angle, A is the azimuth angle, e is the slope of the ground surface, and φ is the tilt azimuth. The pixel value is a digital number proportional to the brightness of one pixel in the image, the pixel value of the Landsat TM image is a value proportional to the luminance of one pixel, and the pixel value of the direct incident illuminance image is This value is proportional to the illuminance of one pixel.
[0012]
The correlation coefficient is an index for evaluating the relationship between the pixel values in the Landsat TM image and the simulation image, and is a numerical value between −1 and +1. The correlation coefficient can be obtained based on a general formula. In order to use the correlation coefficient as an evaluation function, first, a rough shift amount between the Landsat TM image and the simulation image is obtained by visual observation or the like, and one image is centered on the shift amount with respect to the other image. The aforementioned correlation coefficient in a specific pixel region (for example, a region of 100 × 100 pixels at the center of the image) is obtained by moving by a predetermined amount in the pixel direction and the line direction. Then, the amount of movement at which this correlation coefficient is maximized (becomes the optimum value) is determined, the amount of deviation between the two images is obtained from this amount of movement, and the value of the parallel movement parameter is set so that this amount of deviation can be corrected. decide.
[0013]
Directly, p ₀ And l ₀ May be a parameter to be optimized or optimized using an evaluation function. Specifically, however, the above two equations are simplified as follows.
[0014]
p = ax−by + c
l = bx−ay + d
However, a thru | or d are variables. Then, c and d in this equation are set as parallel movement parameters. These parameters are determined so that the correlation coefficient between the pixel values of the two images is an evaluation function so that the evaluation function becomes an optimum value. ₀ And l ₀ Ask for.
[0015]
If geometric transformation is performed for all coordinates of the target area in accordance with the above-described geometric transformation formula using the translation parameter thus determined, high-accuracy geometric correction can be performed over the entire target area without acquiring a large number of GCPs. It can be carried out.
[0016]
If the values of the correlation coefficient in the pixel direction and the correlation coefficient in the line direction when the translation parameter is changed are the maximum values of the translation parameter, fine adjustment within one pixel is performed. Is possible.
[0017]
If the features of the present invention are expressed more abstractly, the features of the present invention are related to the map features obtained in relation to the image features (for example, the aforementioned luminance) appearing in the entire Landsat TM image and the map information on the map. While performing evaluation using an evaluation function that quantitatively evaluates the spatial consistency with the corresponding map information corresponding to the feature on the image (for example, the above-mentioned direct incident illuminance of the sun) of the information without statistical bias Optimize the translation parameter. Here, one of the “evaluation functions for quantitatively evaluating spatial consistency without any bias” is the pixel value based on the brightness on the TM image and the sun on the direct incident illuminance image (simulation image) of the sun. It is a correlation coefficient between the pixel values based on the direct incident illuminance. Other evaluation functions can be used. For example, when the feature on the TM image is the difference between the land and the sea appearing on this image, and the corresponding map information is the difference between the land and the sea shown on the map, Use the identification rate obtained by judging whether the multiple evaluation points randomly selected for each category stratified by distance from the coastline are land or sea based on the corresponding map information. Can do. Note that “stratified” means that the statistical group is divided into several subgroups (layers or categories) based on prior information (distance from the coastline). Even when the identification rate is used as an evaluation function, parameters are determined so that the evaluation function becomes an optimum value or an appropriate value.
[0018]
If the features of the present invention are expressed more abstractly (generally) as a precise geometric correction method for satellite images, it can be expressed as follows.
[0019]
That is, according to the present invention, the coordinates of the satellite image and the coordinates of the map are matched using a geometric transformation formula including a plurality of parameters in order to express the relationship between the coordinates on the map and the coordinates of the satellite image. For the precise geometric correction method of satellite image that corrects the satellite image, it corresponds to the feature on the image among the map features that appear in the entire satellite image and the map related information obtained in relation to the map information on the map. It is characterized in that a plurality of parameters are optimized using an evaluation function that quantitatively evaluates spatial consistency with corresponding map information without statistical deviation.
[0020]
Further, in order to facilitate the parameter determination, the parameter may be determined through two stages. First of all, on the map, whether or not the multiple evaluation points randomly extracted for each category stratified by distance from the coastline on the Landsat TM image (satellite image) is land or sea is determined on the map. Judgment is made based on the land and sea data shown, and the identification rate is obtained. Then, appropriate values of a plurality of parameters are obtained using this identification rate as an evaluation function. Next, this appropriate value is input to the geometric transformation formula, and the satellite image is geometrically corrected to determine the geometrically corrected satellite image. After that, the direct solar irradiance image of the sun created using the digital elevation model is used as a simulation image, and the correlation coefficient of each pixel value of the simulation image and the geometrically corrected satellite image is used as an evaluation function. Perform optimization. In this way, it is not necessary to visually check the deviation between the Landsat TM image (satellite image) and the simulation image, so that parameter determination can be automated.
[0021]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings. FIG. 1 is a diagram for explaining the data processing flow of the method according to the first embodiment of this invention. In this embodiment, Landsat TM images are handled as satellite images. The information on the map projection ancillary record of the reader file attached to the standard processing data included in the system information of the Landsat TM data includes the center coordinates on the map, ie, the UTM coordinates (x ₀ , Y ₀ ), TM image rotation angle (orientation angle θ) to be transmitted with respect to the UTM coordinate system ₀ ), And spatial resolution (Δ). Spatial resolution is the smallest object unit that can be resolved by the sensor that captured the image. Using these pieces of information, the UTM coordinates x and y on an arbitrary map are converted into positions (lines and pixels) on a satellite image (so-called inverse conversion) using the following geometric transformation equations (1) and (2). )can do. These equations (1) and (2) are attached to the system information.
[0022]
p-p ₀ = [Cos θ ₀ * (Xx ₀ ) -Sinθ ₀ * (Y-y ₀ ]] / Δ (1)
l−10 = [− sin θ ₀ * (Xx ₀ ) -Cosθ ₀ * (Y-y ₀ )] / Δ (2)
In the above formula, x ₀ And y ₀ Is the center coordinate of the map coordinate system, x and y are the coordinates of any position in the map coordinate system, and p ₀ And l ₀ Is the center coordinates (pixel number and line number) of the coordinate system of the Landsat TM image, p and l are coordinates (pixel number and line number) of an arbitrary position on the Landsat TM image, and θ ₀ Is the orientation angle and Δ is the spatial resolution.
[0023]
If the above equation is simplified by an affine transformation equation, it becomes as follows.
[0024]
p = ax−by + c
l = bx−ay + d
However, a thru | or d are variables. In this case, c and d are parallel movement parameters. In the system geometric transformation, all the coefficients in the above equation are determined. In the method of the present embodiment, two parameters c and d which mean parallel movement on image coordinates are considered as variables. In addition, p in the above formulas (1) and (2) ₀ And l ₀ It is the same even if is considered as a parameter (variable).
[0025]
In this embodiment, a correlation coefficient between a Landsat TM image that has been subjected to geometric correction through a predetermined geometric transformation (converted satellite image) and a direct incident illuminance image created by simulation is used as an evaluation function for geometric correction. Determine the optimum value of the parameter. Since the correlation coefficient increases as the positioning is correct, it is possible to optimize the parameters of the geometric transformation formula used for the precise geometric correction. Optimization is achieved by searching for parameters so that the correlation coefficient is large.
[0026]
In this embodiment, in order to include the orthographic projection in the geometric transformation using the above formula, the apparent position data on the satellite image obtained by adding the positional deviation based on the undulation to the pixel value obtained by the above formula is used. is doing. An orthographic projection is incorporated in the geometric transformation used to obtain the optimal solution. That is, the positional deviation Δp based on the undulation is added to the pixel value p obtained by substituting the map coordinates (x, y) into the above equations (1) and (2). However, in order to accurately obtain the positional deviation Δp based on the undulations, calculation in consideration of the curvature of the earth is necessary, and the burden on the computer increases. Note that the calculation time can be shortened by utilizing the fact that the positional deviation due to strict calculation is about 10% larger than the positional deviation when the curvature of the earth is not taken into consideration. Further, the positional deviation based on the undulations of the point directly below the satellite may be estimated from the locus of missing data in the original image. In this regard, in “Photogrammetry and Remote Sensing”, Vol. 37, no. 4, page 12-22, published in a paper entitled “Orthographic projection of Landsat TM image and its evaluation”.
[0027]
The direct incident illuminance cos β of the sun necessary for obtaining a simulation image is given by the following equation as a function of the sun position (zenith angle θ and azimuth angle A), the inclination e of the ground surface, and the inclination azimuth φ.
[0028]
cos β = cos θ cos (e) + sin θ sin (e) cos (φ−A)
Here, for the data of the sun position (zenith angle θ and azimuth angle A), the sun position attached to the standard processing data is used. In order to obtain the inclination e and the inclination azimuth φ, a digital elevation model is required. The digital elevation model uses a numerical map 50m (elevation) issued by the Geospatial Information Authority of Japan with projection transformation (UTM coordinate system) and rearrangement (30m). The numerical map 50m mesh (elevation) issued by the Geospatial Information Authority of Japan is composed of lattice points (mesh) composed of equally spaced latitude and longitude, and the size is about 50m. Therefore, in order to use this for the image processing of Landsat TM, projection transformation (UTM coordinate system) and rearrangement (spatial resolution 30 m) are required. In addition to this point in detail Related The Japan Remote Sensing Society Journal, Vol. 21, no. 2, page 150 to 157, is published in a paper entitled “Evaluation Method of Interpolation Method Used for Projection Conversion of Digital Elevation Model”. In the present embodiment, the inclination e and the inclination azimuth φ are obtained from this projection transformation (UTM coordinate system) and rearrangement (spatial resolution 30 m).
[0029]
In the present embodiment, five locations in Iwate Prefecture shown in FIG. 2 (regions in which numbers 1 to 5 are indicated in black circles) are set as target regions. A digital elevation model was created assuming that the size of the satellite image cut out at each point was 430 vertical and 430 horizontal, and the spatial resolution was 30 m.
[0030]
The satellite image is Landsat TM data (pass) taken on June 12, 1995. 107 , Row 32) was used. The solar altitude at the time of photographing was 58 degrees, and the azimuth was 112 degrees clockwise from north. A simulation image of the direct incident illuminance of the sun was created using this information. A simulation image of the point 1 is shown in FIG.
[0031]
First, a rough alignment was performed by visually comparing the satellite image after the system geometric transformation and the simulation image. As a result, a positional shift of 40 pixels in the pixel direction and 27 pixels in the line direction was confirmed. Therefore, the converted Landsat TM image with respect to the simulation image in units of pixels (one pixel at a time) in the pixel direction and the line direction, centered on the positional deviation visually confirmed (40 pixels in the pixel direction and 27 pixels in the line direction). The correlation coefficient between the pixel values of, for example, 100 × 100 pixels located at the center of the Landsat TM image and the pixel values of the corresponding pixels of the simulation image was obtained. The results are shown in Table 1 below.
[0032]
[Table 1]

The numbers on the horizontal axis in Table 1 above mean that the pixels are shifted by ± 2 pixels around 40 pixels in the pixel direction, and the numbers on the vertical axis mean that they are shifted ± 2 pixels around the 27 pixels in the line direction. The numbers in the table are correlation coefficients of pixel values when moving in pixel units. The result of Table 1 shows the evaluation function. From Table 1, it can be seen that the correlation coefficient 0.695 at the position shifted by 40 pixels in the pixel direction and 27 pixels in the line direction is the largest. That is, the optimum value of the evaluation function is to shift 40 pixels in the pixel direction and 27 pixels in the line direction. At this stage, the translation parameters c and d described above can be determined so as to obtain such a deviation amount.
[0033]
FIG. 4 shows a satellite image (band 5) after conversion of point 1. Compared with the simulation image shown in FIG. 3, it can be confirmed that the shadows of both match.
[0034]
Next, a description will be given of a method for determining the optimum value of the parallel movement parameter from the geometric conversion equation with the above-mentioned positional deviation as the origin, which is performed for fine adjustment within one pixel. Pixel values of 100 × 100 pixels in the central region of a satellite image (Landsat TM image shifted by 40 pixels in the pixel direction and 27 pixels in the line direction) created by changing the translation parameters c and d, and the corresponding simulation The correlation coefficient of pixel values of image pixels was calculated. Examples of the results are shown in FIGS. FIG. 5 shows the change in the correlation coefficient when the translation parameter c that determines the movement amount in the pixel direction is changed. The phase when the translation parameter d that determines the movement amount in the line direction is changed is shown in FIG. The change in the number of relationships is shown in FIG. Comparing FIG. 5 and FIG. 6, the change is more obvious when the parameter c is changed than when the parameter d is changed, but both have clean symmetry. By searching two parameters in a two-dimensional manner, an optimum parameter can be determined on the order of an order of magnitude smaller than one pixel. This seems to be because the accuracy of the results is improved because so many points can be used in the calculation of the correlation coefficient. That is, the values of the parameters c and d (0.2 and 0.2 in this example) when the correlation coefficient reaches the maximum value by changing the parameters c and d can be determined as the optimum values. If the image range for calculating the correlation coefficient is a range where the ground cover is uniform and the undulation is relatively large, a correlation coefficient exceeding 0.8 can be obtained. However, since there is no significant change in the positional deviation information even in such a case, it is preferable that the range in which the correlation coefficient is taken in the sense that the number of samples is simply increased.
[0035]
Table 2 shows the result of determining the optimum amount of parallel movement for each of the points 1 to 5 in FIG.
[0036]
[Table 2]

In Table 2, c and d on the horizontal axis indicate parameters, and r indicates the value of the correlation coefficient. Although the selected five points occupy a fairly wide range, the difference between the optimal translation parameters c and d was about 1.1 pixels in both the pixel direction and the line direction. In this range, even if c and d are fixed, they are less than one pixel in the line direction within one pixel. This is considered to indicate that the accuracy of geometric correction based on system information is quite good. In addition, systematic features such as a large shift to the minus side in the pixel direction at points on the west side (points 2 and 5 in FIG. 1) and a large shift to the plus side in the line direction can be seen. This suggests that optimization including rotation is necessary when a larger range is considered.
[0037]
In the above description, the optimum values of the parameters c and d are obtained. However, p in the above formulas (1) and (2) ₀ And l ₀ That is, the center of the scene can be used as a parameter. In this case, the center (p ₀ , L ₀ ) Is changed (δp, δl) to adjust the positional deviation in units smaller than one pixel. In addition, it is expected to produce an output image with high accuracy in a wide area by performing orthographic projection. Furthermore, the optimal solution can be calculated using the Powell method incorporated in IDL (Interactive Data Language) (Ver5.3).
[0038]
In the example of FIG. 1, after the parameters are optimized, the parameters are substituted into the geometric transformation equations (1) and (2) to perform precise geometric transformation of the Landsat TM image. However, as shown in FIG. 7, when transforming Landsat TM data, using equations (1) and (2), the appropriate translation parameter values are first inserted into these equations and It is also possible to execute the optimization and change the translation parameter based on the result.
[0039]
In this embodiment, a simple method (modified system geometric transformation method) for geometric transformation of Landsat TM standard processing data has been proposed. This modified system geometric transformation method considers transformation based on information about map projection included in the system information and taking into account positional displacement (orthogonal projection) due to parallel movement and undulations on the ground. Then, the parallel movement parameter is adjusted using the correlation coefficient between the simulation image of the direct incident illuminance of the sun and the pixel value of the satellite image as an evaluation function. Fine adjustment in units of 0.1 pixel is also possible. In the conventional precise geometric correction using GCP, two-step operation of GCP selection and statistical estimation of the coordinate transformation formula is required, whereas in the method of the present invention, the GCP that has been considered as a problem in the past is required. There is an advantage that selection can be omitted.
[0040]
In the above embodiment, the two translation parameters are optimized. However, the parameters that can be optimized are not limited to the translation parameters, and other parameters are optimized by the same method. Of course, it can be done.
[0041]
FIG. 8 is a diagram illustrating a flow of data processing according to the second embodiment of this invention. In the first embodiment described above, the parameter is optimized after the amount of deviation between the satellite image and the simulation image is visually observed. In this embodiment, instead of detecting the amount of deviation by visual observation, a global amount of deviation is automatically detected. Therefore, in this embodiment, after detecting a global shift amount by a global method (a method using the difference between sea and land), the local method employed in the first embodiment described above is used. The parameter is determined by the method (method using the estimated value of direct incident illuminance of the sun).
[0042]
Specifically, in this embodiment, it is first determined whether a plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline on the Landsat TM image is land or sea. Is determined based on the land and sea data shown on the map, and the identification rate is obtained. Then, an appropriate value of the translation parameter is obtained using this identification rate as an evaluation function. Thereafter, an appropriate value is input to the geometric transformation formula to geometrically correct the Landsat TM image to determine the geometrically corrected Landsat TM image. Next, the solar direct irradiance image created using the digital elevation model is used as a simulation image, and the correlation coefficient of each pixel value of the simulation image and the geometrically corrected Landsat TM image is used as an evaluation function. To do.
[0043]
The global method (method using the difference between sea and land) is described below. The biggest feature that can be identified by satellite images is the difference between land and sea, and these can be clearly identified especially in the near-infrared band images. Whether each point on the map coordinates is sea or land can be determined by numerical information issued by the Geospatial Information Authority of Japan. To use this feature, multiple evaluation points are selected from land and sea. Next, the coordinates on the satellite image of the evaluation point are obtained using the set geometric correction formula. A ratio of evaluation points at which the sea and the land that can be determined from the luminance values on the coordinates coincide with the map information is defined as the fitness. In this method, setting of evaluation points is important. By including not only the location close to the coast but also the location far, optimization can be performed regardless of the initial value of the parameter. In addition, since the calculation time can be shortened by adjusting the evaluation points, it is effective when targeting a wide area such as the entire scene or estimating many parameters.
[0044]
In this embodiment, a satellite image (TM data) “Photo taken on June 12, 1995 (pass 107 / row 32)” from Landsat 5 provided by the Space Development Agency was used. The data used here is so-called standard processing data (level 1b). (1) Sensitivity difference adjustment between sensing elements and (2) Projection to UTM coordinate system using satellite orbit / attitude information And rearrangement (system geometric correction). However, the system geometric correction (1) does not consider the influence of the undulation of the ground surface, (2) the trajectory information is not accurate, (3) the image arrangement does not correspond to the xy axis of the UTM coordinate system ( Due to problems such as orientation angle, it cannot be overlaid with map data accurately.
[0045]
FIG. 9 shows an example of an image file in the visible band. A header record (one line) at the beginning of the file, the first 32 bytes of each record are a record header, and the last 68 bytes are a record trailer, but does not include information necessary for this analysis. The other portions are image data, but the left and right black portions (luminance value 0) are unphotographed portions generated during system correction. In precise geometric transformation (orthogonal projection), an intermediate point of the trajectory of this boundary for each line is set as a point directly below the satellite.
[0046]
In addition to the image file, the standard processing data is accompanied by an information file (lead file). This file contains the following information necessary for geometric transformation.
[0047]
Solar altitude θ = 62 degrees
Solar orientation φ = 119 degrees
Orientation angle θ ₀ = 0.1830542 radian
Scene center (UTM coordinate system)
x ₀ = 534146.1182 m
y ₀ = 4463628.9062 m
Spatial resolution 28.5m
Also, the scene center on this image file is
p ₀ = 3492.0 pixel
l ₀ = 2983.5 line
It is.
[0048]
The position (p, l) on the image when the UTM coordinates (x, y) are given from these pieces of information is the following using the equations (1) and (2) used in the first embodiment. It is calculated by the following formula.
[0049]
p = p ₀ + [Cos θ ₀ * (Xx ₀ ) -Sinθ ₀ * (Y-y ₀ )] / 28.5
l = l ₀ + [-Sinθ ₀ * (Xx ₀ ) -Cosθ ₀ * (Y-y ₀ )] / 28.5
FIG. 10 shows the range of the satellite data on the map. In FIG. 10, a region A is the range of the satellite image of FIG. 9, a region B is a range where the earth method is performed, and a region C is a range where the local method is performed.
[0050]
The data processing flow shown in FIG. 8 can be divided into (A) creation of data for evaluation and (B) estimation of a geometric transformation expression by optimization of the evaluation function. The latter includes a global method using discrimination between sea and land and a local method using direct incident illuminance. Each can be used independently, but in the case of Landsat TM, after creating a geometric transformation formula including rotation by a global method, a local method is used by correcting its translation parameter. Is reasonable.
[0051]
In the data processing flow of FIG. 8, first, data for evaluation is created. There are two types of evaluation data (sea and land discrimination data and sun direct incident illuminance), both of which are created from a numerical map 50m mesh (elevation) published by the Geographical Survey Institute.
[0052]
FIG. 11 shows a cutout of the digital elevation model of a rectangular area B where most of the coastline overlaps with the range of the target satellite image (for path 107 / row 32) shown in FIG. Using this image, a distance image from the coastline (level 1 where the sea and the land touch each other) was created. The area about 750m away from the coastline is divided into 15 levels by distance, and 200 pixels (half sea and land) are randomly extracted for each level (category), and a total of 3000 points are evaluated. Data as shown in Table 3 was created as points. FIG. 12 shows a part of this. In FIG. 12, + marks indicate evaluation points. Further, the position data indicated by the longitude and latitude is converted into the UTM coordinate system as the data used for determining the geometric conversion formula.
[0053]
[Table 3]

FIG. 13 shows a direct incident illuminance image (simulation image) of the sun in the region C shown in FIG. This image is created by projecting / rearranging a numerical elevation model (UTM coordinate system / spatial resolution 30 m) and a numerical map 50 m mesh (elevation) (alternative interpolation method). The size of the image is 430 × 430. The method for obtaining the direct incident illuminance cos β of the sun is the same as that in the first embodiment.
[0054]
Next, estimation of the geometric transformation formula by optimization of the evaluation function will be described.
[0055]
A global method is applied to the band 4 satellite image shown in FIG. The position of the evaluation point on the satellite image is calculated based on the system information, and the case where the luminance value (pixel value) of each evaluation point is equal to or less than the threshold value γ (20 is set) is the sea area, and the luminance value is larger than the threshold value. Judged as land area. And the correctness of this judgment result (identification data) was judged in light of the map relevant information obtained from the above-mentioned numerical map 50m mesh (elevation), and the identification rate was calculated | required. As a result, 718 points out of 2878 points excluding 122 points outside the range of the satellite image from the evaluation points were erroneously determined.
[0056]
In this example, the wrong number is used as an identification rate, and this is used as an evaluation function, so that the evaluation value becomes large (so that the identification rate becomes small), and the optimum center deviation (Δp, Δl) and threshold (γ). As a result, errors were reduced to 51. Here, the optimization is performed by changing Δp, Δl, and the threshold value γ according to a known optimization method, and comparing the brightness value (pixel value) and the threshold value each time. The parameter values at this time are as follows.
[0057]
Case A: Δp = 40.90, Δl = 27.18, γ = 28.39
Case A has only parallel movement, but in order to include rotation, further optimization including the orientation angle deviation Δθ resulted in 40 errors. FIG. 14 shows a wrong place and a point outside the range on the map. Visible in the north are evaluation points that are outside the range of the satellite image. The parameter values at this time are as follows.
[0058]
Case B: Δp = 40.25, Δl = 27.18, Δθ = 0.000217, γ = 29.08
When optimization was performed with the spatial resolution of 28.5 m also changed, no error phenomenon was observed. In addition, although the wrong point is not checked individually, there is no systematic tendency except that many influences of the lakes around Lake Mutsu Ogawara are seen. This may be due to the effects of small islands, estuaries, or port recovery.
[0059]
The images were rearranged using the geometric transformation formula obtained by the global method described above. That is, an approximate amount of deviation between the geometrically corrected Landsat TM image and the simulation image was determined. This rearrangement of images corresponds to the determination of the shift amount of the center of the scene by visual observation in the first embodiment. Thereafter, in the same manner as in the first embodiment, the Landsat TM image geometrically corrected in pixel units is moved in the pixel direction and the line direction, and the phase of each pixel value in the predetermined region of the converted satellite image and simulation image is changed. Change in the number of relationships was sought. Tables 4 and 5 show changes in correlation coefficients obtained when the parameters of cases A and B described above are used, respectively. In this example, in order to avoid an error due to the undulation of the terrain, conversion was performed in consideration of orthographic projection. For interpolation, bilinear interpolation was used. In the analysis, the correlation coefficient obtained by shifting the converted satellite image in the pixel direction and the line direction was used. If the conversion is correct, the maximum value is obtained where there is no deviation.
[0060]
[Table 4]

[Table 5]

In both Table 4 and Table 5, the horizontal axis indicates the amount of movement (number of pixels) in the pixel direction, and the vertical axis indicates the amount of movement (number of pixels) in the line direction. The numbers on the horizontal axis and the vertical axis indicate that the image is shifted in the pixel direction and the line direction in units of pixels on the basis of the deviations Δp and Δl of the center of the scene obtained in Case A and Case B (based on 0) I mean.
[0061]
In both case A (only translation) and case B (including rotation), there is no shift in the line direction, and a maximum shift of about two pixels is seen in the pixel direction. When both are compared, the positional deviation of the latter is slightly smaller.
[0062]
Next, parameters were optimized for each case using the correlation coefficient as an evaluation function. The results are shown below.
[0063]
Case A: Δp = 38.79, Δl = 27.45, correlation coefficient 0.692
Case B: Δp = 39.19, Δl = 27.21, Δθ = −0.000213, correlation coefficient 0.692
When optimization is performed, the effect of including rotation hardly appears. This is because the influence of rotation appears only about a few meters (10000 × Δθ) in the range of the target image (about 10 km). Therefore, a global method is indispensable for obtaining the rotational deviation. In addition, when optimization is performed by a local method, it may be difficult to find a solution unless the initial value is set appropriately. This problem can be solved by using a parameter obtained by a global method as an initial value.
[0064]
【The invention's effect】
According to the present invention, since an appropriate translation parameter can be obtained using the correlation function as an evaluation function, a large number of GCPs can be obtained by performing geometric transformation on all coordinates of the target region using this translation parameter. There is an advantage that the geometric correction can be performed with high accuracy over the entire target area.
[Brief description of the drawings]
FIG. 1 is a diagram for explaining a data processing flow of a method according to a first embodiment of this invention;
FIG. 2 is a diagram showing a target area.
FIG. 3 is a diagram showing a simulation image of point 1 shown in FIG. 2;
4 is a diagram showing a satellite image (band 5) after conversion at the point 1 shown in FIG. 2;
FIG. 5 is a diagram illustrating a change in correlation coefficient with respect to a parameter c.
FIG. 6 is a diagram illustrating a change in correlation coefficient with respect to a parameter d.
FIG. 7 is a diagram illustrating another data processing flow according to the first embodiment;
FIG. 8 is a diagram illustrating a flow of data processing according to the second embodiment of this invention;
FIG. 9 is a diagram illustrating an image of an image file in a visible band.
FIG. 10 is a diagram showing a target area.
11 is a diagram showing a digital elevation model of a region B in FIG.
FIG. 12 is a diagram showing an example of evaluation points randomly extracted.
13 is a diagram showing a direct incident illuminance image (simulation image) in region C shown in FIG.
FIG. 14 is a diagram showing evaluation points with erroneous land and sea determinations.
[Explanation of symbols]
A to C region

Claims (4)

In order to match the coordinates of the Landsat TM image with the coordinates of the map using a geometric transformation formula including a plurality of parameters in order to express the relationship between the coordinates on the map and the coordinates of the Landsat TM image. A precision geometric correction method for a Landsat TM image for correcting a Landsat TM image,
Using a geometric transformation formula including a translation parameter included in system information provided accompanying the Landsat TM image as the geometric transformation formula,
The spatial consistency between the feature on the image appearing in the entire Landsat TM image and the corresponding map information corresponding to the feature on the image among the map related information obtained in association with the map information on the map is statistically calculated. using an evaluation function to evaluate quantitatively without deviation, we have rows optimization of the translation parameters,
The feature on the image is the difference between the land and the sea appearing on the image,
The correspondence map information is a difference between the land and the sea shown on the map,
As the evaluation function, whether or not the identification data indicating whether the plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline is land or sea is determined based on the corresponding map information. precise geometric correction method features and to Lula Ndosatto TM images using the identification rate. The geometric transformation formula included in the system information is:
p−p ₀ = [cos θ ₀ * (x−x ₀ ) −sin θ ₀ * (y−y ₀ )] / Δ
l−l ₀ = [− sin θ ₀ * (x−x ₀ ) −cos θ ₀ * (y−y ₀ )] / Δ
Where x ₀ and y ₀ are center coordinates of the coordinate system of the map, x and y are coordinates of an arbitrary position of the coordinate system of the map, and p ₀ and l ₀ are the Landsat The center coordinates of the TM image coordinate system, p and l are the coordinates of an arbitrary position on the Landsat TM image, θ ₀ is the rotation angle of the image coordinate system with respect to the map coordinate system, and Δ Is the spatial resolution of the Landsat TM image,
The above equation, p = ax-by + c and l = bx-ay + d and when simplified (but a to d are variables), the c and d in claim 1, characterized in that said translation parameters A precision geometric correction method for the described Landsat TM image. The geometric transformation formula included in the system information is:
p−p ₀ = [cos θ ₀ * (x−x ₀ ) −sin θ ₀ * (y−y ₀ )] / Δ
l−l ₀ = [− sin θ ₀ * (x−x ₀ ) −cos θ ₀ * (y−y ₀ )] / Δ
Where x ₀ and y ₀ are center coordinates of the coordinate system of the map, x and y are coordinates of an arbitrary position of the coordinate system of the map, and p ₀ and l ₀ are the Landsat The center coordinates of the TM image coordinate system, p and l are the coordinates of an arbitrary position on the Landsat TM image, θ ₀ is the rotation angle of the image coordinate system with respect to the map coordinate system, and Δ Is the spatial resolution of the Landsat TM image,
Precise geometric correction method of Landsat TM images according to claim 1, wherein said p ₀ and l ₀ is a parameter for optimizing using the evaluation function. The satellite image is used to match the coordinates of the satellite image with the coordinates of the map using a geometric transformation formula including a plurality of parameters to express the relationship between the coordinates on the map and the coordinates of the satellite image. A precise geometric correction method for a satellite image for geometric correction,
A plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline on the satellite image indicates whether land or sea is judged whether it is land or sea. Judgment based on the data, to determine the identification rate,
Obtaining appropriate values of the plurality of parameters using the identification rate as an evaluation function
Input the appropriate value to the geometric transformation formula to geometrically correct the satellite image to determine a geometrically corrected satellite image;
The solar direct irradiance image created using a digital elevation model is used as a simulation image.
An accurate geometric correction method for a satellite image, wherein the plurality of parameters are optimized by using a correlation coefficient between pixel values of the simulation image and the geometrically corrected satellite image as an evaluation function.

Images