画素値補正プログラムおよび記録媒体
【課題】近赤外バンドを含む光学センサにより取得した起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くならず大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去できる物理的に信頼性の高い画素値補正プログラム等を提供する。
【解決手段】近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する。同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する。数値標高モデルの座標データから分類クラス毎に各画素の標高を抽出し、分類クラス毎に分析ステップにより画素毎に大気影響度を求める。各画素の大気影響度に対して中央値フィルタを適用し画素毎に地域的な大気影響度を求め、各画素の画素値及び標高と基準標高と地域的な大気影響度とに基づき各画素の画素値を補正する。
【解決手段】近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する。同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する。数値標高モデルの座標データから分類クラス毎に各画素の標高を抽出し、分類クラス毎に分析ステップにより画素毎に大気影響度を求める。各画素の大気影響度に対して中央値フィルタを適用し画素毎に地域的な大気影響度を求め、各画素の画素値及び標高と基準標高と地域的な大気影響度とに基づき各画素の画素値を補正する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムおよび記録媒体に関する。
【背景技術】
【０００２】
地球観測衛星は種々のセンサを搭載しており、収集された衛星画像データは気象現象、陸域の植生被覆、水資源の調査等様々な目的に利用されている。例えば、地球観測衛星の１つであるランドサット（Landsat）４号には７バンドのＴＭ(Thematic Mapper)センサが搭載されており、バンド１〜３は可視バンド、バンド６は温度情報を得るための熱バンド、バンド４、５および７は近赤外バンドとなっている。
【０００３】
衛星画像データの輝度値（画素値）には、地表で反射された光等以外に、大気によって散乱された散乱光であるパスラディアンス（path radiance : 光路輝度）が影響を与えている。このパスラディアンスの特性は、大気中のもや、霞または埃のヘイズ（haze）または霧を含む大気状態により異なる。可視バンドから中間赤外バンドのように近赤外バンドを含む光学センサによって取得された衛星画像データには、大気と地形との双方の影響に基づく変動が含まれている。特に、起伏の大きい地域の可視バンド衛星画像データには大気状態の空間的変動に起因するパスラディアンスが大きく影響しているため、この影響を除去することは重要である。
【０００４】
非特許文献１では、近赤外バンドにおける各画素は可視バンドにおけるある画素に対応するものと仮定し、散乱が可視バンドで生じている場合は近赤外バンドにおける各画素は可視バンドにおける複数の画素に対応しているものと考えている。そこで、近赤外バンドで土地被覆分類等による分類により同じ分類クラスに分類された可視バンドの各画素の画素値を、それらの画素値の平均値で置き換えるという方法を大気補正法として提案している。この大気補正法は局所的にヘイズがある衛星画像データに対しては見かけ上良くなるが、局所的な情報は失われることになるという問題があった。例えば、近赤外バンドではヘイズのある地域での水田もクリアな地域での水田も同じ分類クラスに属しているが、可視バンドでは両者は異なる分類クラスに属している。ここで水田がヘイズのある地域に偏在していると、可視バンドにおける画素値の平均値はヘイズのある地域側に偏った値となる。このため、クリアな地域の水田の画素値を当該平均値で置き換えると、クリアな地域の水田も汚されてしまうことになる。つまり、非特許文献１の大気補正法では水田がヘイズのある地域に偏在しているというような局所的な情報は失われてしまうことになるという問題があった。さらに、非特許文献１では起伏のある地域を扱っていないため、起伏による陰影の情報という局所的な情報も失われてしまい、見かけも良くなくなってしまうという問題があった。
【０００５】
非特許文献２では、起伏の激しい山岳地域では地表面の傾斜による太陽入射照度の違いが重要な問題となる点に着目し、大気・地形効果補正の簡便な方法を提案している。詳しくは、パスラディアンスが標高の一次関数として表すことができるということを放射伝達のシミュレーションと実際の衛星画像データとを用いて示している。しかし、大気状態が水平方向に変化する場合は取り扱っていないという問題があった。
【０００６】
非特許文献３では、タセルドキャップ（Tasseled Cap）というランドサットＴＭの上述の各バンドの画素値を線形変換した指標の中の４番目の指標（簡便にはバンド１とバンド３との差）がヘイズの量と相関が高いという経験式を用いた方法が記載されている。しかし、非特許文献３では画素毎に分光反射特性が異なるということが無視されている。統計的にも変換後のバンド間の相関が高くなることが指摘されており、非特許文献３の方法は物理的に信頼性が低いという問題があった。
【０００７】
非特許文献４では、まず赤外バンドを分類し、その分類クラス毎にエアロゾル（ヘイズ）の光学的厚さを求め、次に、当該光学的厚さを空間的に平均した値を求めている。詳しくは、光学的厚さの空間的な平滑化を比較的狭い範囲（１５０ｍ程度）で通常の平均値を用いて行なっている。この平均値を利用して、画素値から地表面の反射率を補正している。しかし、非特許文献４では起伏のある地域を取り扱っていないという問題があった。さらに、光学的厚さおよび反射率という物理量を用いているため、クリアな地域における大気の散乱等を予め計算しておく必要がある。従って、クリアな地域を厳密に定めておく必要があり、煩雑で非実用的であるという問題があった。加えて、補正に利用する値は比較的狭い範囲における通常の平均値であるため、大きな範囲で平滑化を行った場合、必要とする変化を捉えることができないという問題があった。
【０００８】
【非特許文献１】M.J.Carlooto:“Reducingthe effects of space-varying, wavelength-dependent scattering in multispectralimagery”, Int.J.Remote Sensing, 1999, Vol.20, No.17, pp.3333-3344.
【非特許文献２】飯倉善和、横山隆三、「ランドサットＴＭ画像の大気および地形効果の補正」、日本リモートセンシング学会誌、１９９９年、第１９巻、第１号、ｐ．２−１６．(Y.Iikura and R.Yokoyama:“Correction of Atmospheric and TopographicEffects on Landsat TM Imagery”, Journal of Remote Sensing Society of Japan, 1999,Vol.19, No.1, pp.2-16.)
【非特許文献３】J.Lavreau:“De-hazingLandsat Thematic Mapper images”, Photogrammetric Engineering & RemoteSensing, 1991, Vol.57, No.10, pp.1297-1302.
【非特許文献４】S.Liang,H.Fang, and M.Chen:“Atmospheric Correction of Landsat ETM+ Land Surface Imagery- Part I：Methods”, IEEE Trans. on Geoscience and Remote Sensing, 2001, Vol.39,No.11, pp.2490-2498.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
上述のように、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データには大気状態の空間的変動に起因するパスラディアンスが大きく影響しているため、この影響を除去することは重要である。しかし、非特許文献１に記載された大気補正法では、局所的にヘイズがある衛星画像データに対しては見かけ上良くなるが、局所的な情報は失われることになるという問題があった。さらに、非特許文献１では起伏のある地域を扱っていないため、起伏による陰影の情報という局所的な情報も失われてしまい、見かけも良くなくなってしまうという問題があった。非特許文献２に記載された大気・地形効果補正の方法では、起伏のある地域を考慮しているが、大気状態が水平方向に変化する場合は取り扱っていないという問題があった。非特許文献３の方法は物理的に信頼性が低いという問題があった。非特許文献４に記載された方法では、起伏のある地域を取り扱っていないという問題があった。さらに、光学的厚さおよび反射率という物理量を用いているため、クリアな地域における大気の散乱等を予め計算しておく必要がある。従って、クリアな地域を厳密に定めておく必要があり、煩雑で非実用的であるという問題があった。加えて、補正に利用する値は比較的狭い範囲における通常の平均値であるため、大きな範囲で平滑化を行った場合、必要とする変化を捉えることができないという問題があった。
【００１０】
そこで、本発明の目的は、上記問題を解決するためになされたものであり、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することにある。さらに、本発明の目的は、簡便で実用的であり、大きな範囲で平滑化を行った場合でも、必要とする変化を捉えることができる画素値補正プログラム等を提供することにある。
【課題を解決するための手段】
【００１１】
この発明の画素値補正プログラムは、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムであって、コンピュータに、前記対象地域に関する近赤外バンドの衛星画像データを入力し、該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ、前記対象地域に関する可視バンドの衛星画像データを入力し、前記分類ステップで同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ、前記対象地域に関する数値標高モデルの座標データを入力し、前記分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ、前記分類クラス毎に、前記可視バンドデータ抽出ステップで抽出された画素値、前記数値標高データ抽出ステップで抽出された標高及び所定の基準標高に基づき画素毎に大気影響度を求める分析ステップ、前記分析ステップで求めた各画素の大気影響度に対して中央値フィルタを適用し、該画素毎に地域的な大気影響度を求めるフィルタステップ、前記可視バンドデータ抽出ステップで抽出した各画素の画素値と、前記数値標高データ抽出ステップで抽出した標高と、前記所定の基準標高と、前記フィルタステップで求めた該画素の地域的な大気影響度とに基づき、各画素の画素値を補正する補正ステップを実行させるための画素値補正プログラムである。
【００１２】
ここで、この発明の画素値補正プログラムにおいて、前記所定の基準は、近赤外バンドの衛星画像データの画素値が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準とすることができる。
【００１３】
ここで、この発明の画素値補正プログラムにおいて、前記所定の基準標高は前記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値とすることができる。
【００１４】
ここで、この発明の画素値補正プログラムにおいて、前記分析ステップは、前記分類クラス毎に、横軸を標高、縦軸を画素値とするグラフを用い各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットするグラフ作成ステップと、クリアな領域における傾きＢ０（画素値／標高）を設定する傾き設定ステップと、前記傾き設定ステップで設定された傾きＢ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算する仮の画素値計算ステップと、前記仮の画素値計算ステップで計算された仮の画素値Ｙｐｉの中央値を求め該中央値を基準標高Ｘ０における画素値Ｙ０とする画素値決定ステップと、前記画素値決定ステップで決定された画素値Ｙ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算し該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とするステップとを備えることができる。
【００１５】
ここで、この発明の画素値補正プログラムにおいて、前記フィルタステップは、各画素の周囲６００ｍ程度以内の画素の大気影響度を抽出し、該大気影響度の中央値を用いて平滑化した結果を各画素の地域的な大気影響度として求めることができる。
【００１６】
ここで、この発明の記載の画素値補正プログラムにおいて、前記補正ステップは、各画素の標高Ｘｉ、画素値Ｙｉ、基準標高Ｘ０及び地域的な大気影響度Ｂｉ’を用いて補正後の画素値Ｙｉ’を、Ｙｉ’＝Ｙｉ―Ｂｉ’×（Ｘ０−Ｘｉ）、ｉ＝１〜全画素数、により計算することができる。
【００１７】
この発明の記録媒体は、本発明の画素値補正プログラムを記録したコンピュータ読取り可能な記録媒体である。
【発明の効果】
【００１８】
本発明の画素値補正プログラムによれば、大気影響度を土地被覆分類等の分類クラス毎に標高を加味した式により計算しているため、タセルドキャップよりも物理的に信頼性の高い画素値補正プログラム等を提供することができる。画素毎の大気影響度から地域的な大気影響度を求める際に中央値フィルタを利用しているため、ヘイズ等の非連続な大気の状況の変化を自然に把握することができ、大きい範囲で平滑化を行った場合であっても、必要な変化を捉えることができる。さらに、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めた地域的な大気影響度を求めることができる。画素における大気影響度または地域的な大気影響度を基準標高における画素に対する傾きと考えたことにより、補正量を標高の一次式として簡単に求めることができる。
【００１９】
この結果、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することができるという効果がある。さらに、光学的厚さまたは反射率等の物理量を用いていないため、簡便で実用的な画素値補正プログラム等を提供することができるという効果がある。
【発明を実施するための最良の形態】
【００２０】
まず、本発明の、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムの基盤となる理論について説明し、次に各実施例について図面を参照して詳細に説明する。
【００２１】
理論の説明．
１．背景技術の非特許文献２で説明した、パスラディアンスを標高の一次関数として表すことができるという点を拡張して、当該一次関数の傾きが場所により異なるという以下の式1および２を本発明の画素値補正プログラムの基本モデルとした。これにより、大気の補正を見通しのよいものとした。標高（ｍ）をＸ軸、画素値（Digital Number :ＤＮ）をＹ軸とし、標高Ｘｉにおける補正前の画素値をＹｉ、補正量をＤｉ、傾きをＢｉ、パラメータをＡ、補正された画素値をＹｉ’とすると、
【００２２】
Ｙｉ’＝Ｙｉ−Ｄｉ （１）
【００２３】
と表すことができる。補正量Ｄｉは、
【００２４】
Ｄｉ＝Ａ−Ｂｉ×Ｘｉ （２）
【００２５】
である。ここで、大気の影響が無視できる標高（基準標高Ｘ０）での補正量を０とすれば、Ａは以下のように設定できる。
【００２６】
０＝Ａ−Ｂｉ×Ｘ０
∴Ａ＝Ｂｉ×Ｘ０
【００２７】
上記パラメータＡを式１に代入すると、式１は、以下の式３のようになり、パラメータはＢｉのみでよいことになる。
【００２８】
Ｄｉ＝Ｂｉ×Ｘ０−Ｂｉ×Ｘｉ＝Ｂｉ×（Ｘ０−Ｘｉ） （３）
【００２９】
対象画素（Ｘｉ，Ｙｉ）の傾きＢｉは、基準標高Ｘ０における画素値Ｙ０（大気の影響がない場合の画素値）がわかれば、以下の式４から求めることができる。
【００３０】
Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ） （４）
【００３１】
２．基本モデルのＢｉを画素（Ｘｉ，Ｙｉ）毎に求めるためには、式４の計算を行うことになる。ここで、式４における画素値Ｙ０は土地被覆分類等により分類された分類クラス毎に必要である。本発明の画素値補正プログラムでは、画素値Ｙ０を求めるために、背景技術の非特許文献１で説明したように近赤外バンドで同じ分類クラスに分類された可視バンドの各画素の画素値を利用する。非特許文献１では可視バンドの各画素の平均値に各画素値を置き換えているが、平均値からの差が大気の影響度と考えることもできる。画像の半分以上がクリアな領域である場合、平均値ではなくそれらの中央値に画素値を置き換えた方が良い。画素値Ｙ０を求めるためにクリアな領域での傾きＢ０を利用する。Ｂ０は６Ｓ（Second Simulation of the Satellite Signal in the Solar Spectrumの略で、フランスのLille科学技術大学で開発されたプログラム。）等により設定した値（バンド１で７〜８ＤＮ／ｋｍ）を用いる。まず、同じ分類クラスに含まれるすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを以下の式５を用いて計算する。
【００３２】
Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ） （５）
【００３３】
クリアな領域に含まれる画素のＹｐｉは標高によらず一定となるため、これを大気の影響がない場合における画素値Ｙ０とみなすことができる。分類クラスの中でクリアな領域に含まれる画素が多い場合、Ｙｐｉの中央値はＹ０となる。クリアでない領域に含まれる画素のＹｐｉは大きくなるため、このような画素が多い場合には、すべての画素（Ｘｉ，Ｙｉ）について計算されたＹｐｉの中央値よりもクリアな地域について計算された中央値を使用すればよい。具体的には、クリアな地域が占める画像中の割合をＷ（％）とした場合、すべての画素（Ｘｉ，Ｙｉ）について計算されたＹｐｉの度数分布の下からＷ／２（％）のＹｐｉを用いて計算された中央値を使用すればよい。
【００３４】
３．上述した１および２により画素毎に傾きＢｉを求めることができる。傾きＢｉは大気の影響度として採用することができる。但し、傾きＢｉには種々の誤差要因が含まれている。衛星画像データからは、大気の影響度には地域的な特徴があり、ある場所ではクリアであるが他の場所では薄い霧またはヘイズがかかっていることがわかる。大気の影響度はある程度の広がりで緩やかに変化すると考えられる。従って、大気の影響度としては画素単位ではなくその近傍のデータも利用した方が誤差を少なくすることができると考えられる。そこで、対象画素の周囲の画素の傾きＢｉを抽出してそれらの中央値を利用するという空間的な中央値フィルタを適用して、地域的な大気影響度Ｂｉ’を求める。平均値ではなく中央値を利用する理由は、クリアな領域とクリアでない領域とがはっきりと分かれている場合、中央値に比べて平均値を利用すると両領域の境界がぼけてしまうためである。この地域的な大気影響度Ｂｉ’を用いて以下の式６により補正量Ｄｉを求める。
【００３５】
Ｄｉ＝Ｂｉ’×（Ｘ０−Ｘｉ） （６）
【００３６】
補正された画素値Ｙｉ’は、式１および６により、以下の式７のようになる。
【００３７】
Ｙｉ’＝Ｙｉ−Ｂｉ’×（Ｘ０−Ｘｉ） （７）
【実施例１】
【００３８】
図１は、本発明の実施例１における画素値補正プログラムの流れをフローチャートで示す。以下では図１における各処理ブロックの説明を、図２〜４の実例、図５に示される細部のフローチャート、図６に示されるグラフ、図７〜９の実例を用いて説明する。
【００３９】
図１に示されるように、コンピュータに、まず、所望の対象地域（例えば岩手県大川地区）に関する近赤外バンドの衛星画像データを入力し、当該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ（ステップＳ１０）を実行させる。入力する近赤外バンドの衛星画像データは一般的なものであるため、説明は省略する。所定の基準としては、例えば近赤外バンドの衛星画像データの画素値（ＤＮ）が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準を用いることができる。他の基準を用いてもよいことは勿論である。
【００４０】
図２は、分類ステップ（ステップＳ１０）を実行させた結果の分類データ画像を示す。図２に示される分類データ画像はカラー合成画像である。
【００４１】
次に、図１に示されるように、コンピュータに、上記対象地域と同じ対象地域に関する可視バンドの衛星画像データを入力し、分類ステップ（ステップＳ１０）で同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ（ステップＳ２０）を実行させる。
【００４２】
図３は、可視バンド分類ステップ（ステップＳ２０）で入力する可視バンド（バンド１）の衛星画像データの一例を示す。図３の右上部Ｈａｚｅで示されるように、当該部分はヘイズの影響により明るくなっている。
【００４３】
続いて、図１に示されるように、コンピュータに、上記対象地域と同じ対象地域に関する数値標高モデルの座標データを入力し、分類ステップ（ステップＳ１０）で分類された分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ（ステップＳ３０）を実行させる。
【００４４】
図４は、画素値取得ステップ（ステップＳ３０）で入力する数値標高モデルの座標データ画像の一例を示す。当該座標データ画像は、国土地理院の数値データ５０ｍメッシュ（標高）からＵＴＭ座標系（Universal Transverse Mercator system）へ変換して作成したものである。
【００４５】
次に、図１に示されるように、コンピュータに、上記分類クラス毎に、可視バンドデータ抽出ステップ（ステップＳ２０）で抽出された画素値と、数値標高データ抽出ステップ（ステップＳ３０）で抽出された標高と、所定の基準標高とに基づき画素毎に大気影響度を求める分析ステップ（ステップＳ４０）を実行させる。
【００４６】
図５は分析ステップ（ステップＳ４０）の詳細な流れをフローチャートで示し、図６は分析ステップ（ステップＳ４０）の詳細を説明するためのグラフを示す。以下の処理は上記分類クラス毎に実行される。図６に示される横軸を標高（ｍ）、縦軸を画素値（ＤＮ）とするグラフを用い、各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットする（グラフ作成ステップ。ステップＳ４１）。
【００４７】
次に、クリアな領域における傾きＢ０（画素値／標高）を設定する（傾き設定ステップ。ステップＳ４２）。上述のように、Ｂ０は６Ｓ等の放射伝達シミュレーションにより設定した値（バンド１で７〜８ＤＮ／ｋｍ）を用いればよい。
【００４８】
傾き設定ステップ（ステップＳ４２）で設定された傾きＢ０を用い、グラフ作成ステップ（ステップＳ４１）でプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、上述の式５、
【００４９】
Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜上記分類クラスの画素数 （５）
【００５０】
を用いて計算する（仮の画素値計算ステップ。ステップＳ４３）。所定の基準標高Ｘ０としては、上記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値とすることが好適である。図６では所定の基準標高Ｘ０として１２００ｍを用いているが、所定の基準標高Ｘ０は１２００ｍに限定されるものではない。
【００５１】
仮の画素値計算ステップ（ステップＳ４３）で計算された仮の画素値Ｙｐｉの中央値を求め、当該中央値を基準標高Ｘ０における画素値Ｙ０とする（画素値決定ステップ。ステップＳ４４）。図６では、画素値Ｙ０は約７５（ＤＮ）となっている。
【００５２】
画素値決定ステップ（ステップＳ４４）で決定された画素値Ｙ０を用い、グラフ作成ステップ（ステップＳ４１）でプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、上述の式４、
【００５３】
Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜上記分類クラスの画素数 （４）
【００５４】
を用いて計算し、当該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とする（ステップＳ４５）。図６では画素（４５０、９０）について基準標高１２００ｍにおける画素（１２００、７５）に対する傾きＢｉが直線ｉにより示されている。この傾きＢｉ、即ち７５０ｍ（＝１２００ｍ−４５０ｍ）当たり１５ＤＮ（＝９０ＤＮ−７５ＤＮ）変わるという値が画素（４５０、９０）における大気影響度である。図７は、上述のようにして求められた画素毎の大気影響度Ｂｉを示す画像である。
【００５５】
図１に戻り、続いてコンピュータに、分析ステップ（ステップＳ４０）で求めた各画素（Ｘｉ，Ｙｉ）の大気影響度Ｂｉに対して中央値フィルタを適用し、画素（Ｘｉ，Ｙｉ）毎に地域的な大気影響度Ｂｉ’を求めるフィルタステップ（ステップＳ５０）を実行させる。図８は、フィルタステップ（ステップＳ５０）の実行結果を示す画像である。図８では、フィルタステップ（ステップＳ５０）において、各画素（Ｘｉ，Ｙｉ）の周囲６００ｍ程度以内の画素の大気影響度Ｂｉを抽出し、この大気影響度Ｂｉの中央値を用いて平滑化した結果を各画素（Ｘｉ，Ｙｉ）の地域的な大気影響度Ｂｉ’として求めている。
【００５６】
最後に、コンピュータに、可視バンドデータ抽出ステップ（ステップ２０）で抽出した各画素の画素値Ｙｉと、数値標高データ抽出ステップ（ステップＳ３０）で抽出した標高Ｘｉと、所定の基準標高Ｘ０と、フィルタステップ（ステップ５０）で求めた当該画素の地域的な大気影響度Ｂｉ’とに基づき、上述の式６、
【００５７】
Ｄｉ＝Ｂｉ’×（Ｘ０−Ｘｉ） （６）
【００５８】
により補正量Ｄｉを求め、上述の式７、
【００５９】
Ｙｉ’＝Ｙｉ−Ｂｉ’×（Ｘ０−Ｘｉ） （７）
【００６０】
により各画素の画素値Ｙｉを補正してＹｉ’を求める補正ステップ（ステップ６０）を実行させる。ここで、ｉ＝１〜全画素数である。図９は、上述した補正ステップ（ステップＳ６０）を実行した結果を示す画像である。
【００６１】
以上より、本発明の実施例１によれば、コンピュータに、まず、所望の対象地域に関する近赤外バンドの衛星画像データを入力し、当該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップを実行させる。次に、上記対象地域と同じ対象地域に関する可視バンドの衛星画像データを入力し、分類ステップ（ステップＳ１０）で同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ（ステップＳ２０）を実行させる。続いて、上記対象地域と同じ対象地域に関する数値標高モデルの座標データを入力し、上記分類クラス毎に、各画素の標高Ｘｉを抽出する数値標高データ抽出ステップ（ステップＳ３０）を実行させる。上記分類クラス毎に、可視バンドデータ抽出ステップ（ステップＳ３０）で抽出された画素値と、数値標高データ抽出ステップ（ステップＳ３０）で抽出された標高と、所定の基準標高とに基づき画素毎に大気影響度を求める分析ステップ（ステップＳ４０）を実行させる。分析ステップ（ステップＳ４０）で求めた各画素（Ｘｉ，Ｙｉ）の大気影響度Ｂｉに対して中央値フィルタを適用し、画素（Ｘｉ，Ｙｉ）毎に地域的な大気影響度Ｂｉ’を求めるフィルタステップ（ステップＳ５０）を実行させる。可視バンドデータ抽出ステップ（ステップＳ２０）で抽出した各画素の画素値Ｙｉと、数値標高データ抽出ステップ（ステップＳ３０）で抽出した標高と、所定の基準標高Ｘ０と、フィルタステップ（ステップ５０）で求めた当該画素の地域的な大気影響度Ｂｉ’とに基づき、上述の式６により補正量Ｄｉを求め、上述の式７により各がその画素値Ｙｉを補正してＹｉ’を求める補正ステップ（ステップＳ６０）を実行させる。
【００６２】
本発明の画素値補正プログラムでは、大気影響度Ｂｉを土地被覆分類等の分類クラス毎に標高Ｘｉを加味した式４により計算しているため、タセルドキャップよりも物理的に信頼性の高い画素値補正プログラム等を提供することができる。本発明の画素値補正プログラムでは、画素毎の大気影響度Ｂｉから地域的な大気影響度Ｂｉ’を求める際に中央値フィルタを利用しているため、ヘイズ等の非連続な大気の状況の変化を自然に把握することができ、大きい範囲で平滑化を行った場合であっても、必要な変化を捉えることができる。さらに、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めた地域的な大気影響度Ｂｉ’を求めることができる。本発明の画素値補正プログラムでは、画素（Ｘｉ，Ｙｉ）における大気影響度Ｂｉまたは地域的な大気影響度Ｂｉ’を基準標高Ｘ０における画素（Ｘ０，Ｙ０）に対する傾きＢｉと考えたことにより、補正量Ｄｉを標高Ｘｉの一次式として簡単に求めることができる。以上より、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することができる。さらに、光学的厚さまたは反射率等の物理量を用いていないため、簡便で実用的な画素値補正プログラム等を提供することができる。
【実施例２】
【００６３】
図１０は、上述した実施例１を実現するための本発明のコンピュータ・プログラム（画素値補正プログラム）を実行するコンピュータの内部回路１０を示すブロック図である。図１０に示されるように、ＣＰＵ１１、ＲＯＭ１２、ＲＡＭ１３、画像制御部１４、コントローラ１８、入力制御部２０および外部インタフェース（Ｉ／Ｆ）部２１はバス２２に接続されている。図１０において、上述の本発明の画素値補正プログラムは、ＲＯＭ１２、ＨＤＤ１７ａまたはＤＶＤ１７ｎ（あるいはＣＤ−ＲＯＭ）等の記録媒体（脱着可能な記録媒体を含む）に記録されている。ＨＤＤ１７ａまたはＤＶＤ１７ｎ等には入力した近赤外バンドの衛星画像データ、可視バンドの衛星画像データ、数値標高モデルの座標データ等が記録されている。本発明の画素値補正プログラムは、ＲＯＭ１２からバス２２を介し、あるいはＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体からコントローラ１８を経由してバス２２を介しＲＡＭ１３へロードされる。入力制御部２０はマウス、テンキー等の入力操作部１９と接続され入力制御等を行う。画像メモリであるＶＲＡＭ１５は、分類データ画像、画素毎の大気影響度Ｂｉを示す画像、フィルタステップ（ステップＳ５０）の実行結果を示す画像、補正ステップ（ステップＳ６０）を実行した結果を示す画像等を表示するディスプレイ等の表示部１６の少なくとも一画面分のデータ容量に相当する容量を有しており、画像制御部１６はＶＲＡＭ１５のデータを画像データへ変換して表示部１６へ送出する機能を有している。外部Ｉ／Ｆ部２１は、近赤外バンドの衛星画像データ、可視バンドの衛星画像データ、数値標高モデルの座標データ等を入力する際におけるインタフェース機能を有する。
【００６４】
上述のようにＣＰＵ１１が上述の本発明の画素値補正プログラムを実行することにより、本発明の目的を達成することができる。当該画素値補正プログラムは上述のようにＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体の形態でコンピュータＣＰＵ１１に供給することができ、当該画素値補正プログラムを記録したＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体も同様に本発明を構成することになる。当該画素値補正プログラムを記録した記録媒体としては上述された記録媒体の他に、例えばメモリ・カード、メモリ・スティック等を用いることができる。
【産業上の利用可能性】
【００６５】
本発明の活用例として、ランドサットＴＭセンサにより取得された衛星画像データの画素値の補正に適用することができる。
【図面の簡単な説明】
【００６６】
【図１】本発明の実施例１における補正プログラムの流れを示すフローチャートである。
【図２】分類ステップ（ステップＳ１０）を実行させた結果の分類データ画像を示す図である。
【図３】可視バンド分類ステップ（ステップＳ２０）で入力する可視バンド（バンド１）の衛星画像データの一例を示す図である。
【図４】分析ステップ（ステップＳ４０）で入力する数値標高モデルの座標データ画像の一例を示す図である。
【図５】分析ステップ（ステップＳ４０）の詳細な流れを示すフローチャートである。
【図６】分析ステップ（ステップＳ４０）の詳細を説明するためのグラフである。
【図７】分析ステップ（ステップＳ４０）で求められた画素毎の大気影響度Ｂｉを示す画像である。
【図８】フィルタステップ（ステップＳ５０）の実行結果を示す画像である。
【図９】補正ステップ（ステップＳ６０）を実行した結果を示す画像である。
【図１０】本発明の実施例を実現するための本発明の画素値補正プログラムを実行するコンピュータの内部回路１０を示すブロック図である。
【符号の説明】
【００６７】
Ｂｉ、Ｂ０ 傾き、 Ｈａｚｅ ヘイズによる影響がある領域、 Ｘ０ 基準標高、 Ｙ０ 基準標高における画素値、 １０ 内部回路、 １１ ＣＰＵ、 １２ ＲＯＭ、 １３ ＲＡＭ、 １４ 画像制御部、 １５ ＶＲＡＭ、 １６ 表示部、 １７ａ ＨＤＤ、 １７ｎ ＤＶＤ、 １８ コントローラ、 １９ 入力操作部、 ２０ 入力制御部、 ２１ 外部Ｉ／Ｆ部、 ２２ バス。
【技術分野】
【０００１】
本発明は、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムおよび記録媒体に関する。
【背景技術】
【０００２】
地球観測衛星は種々のセンサを搭載しており、収集された衛星画像データは気象現象、陸域の植生被覆、水資源の調査等様々な目的に利用されている。例えば、地球観測衛星の１つであるランドサット（Landsat）４号には７バンドのＴＭ(Thematic Mapper)センサが搭載されており、バンド１〜３は可視バンド、バンド６は温度情報を得るための熱バンド、バンド４、５および７は近赤外バンドとなっている。
【０００３】
衛星画像データの輝度値（画素値）には、地表で反射された光等以外に、大気によって散乱された散乱光であるパスラディアンス（path radiance : 光路輝度）が影響を与えている。このパスラディアンスの特性は、大気中のもや、霞または埃のヘイズ（haze）または霧を含む大気状態により異なる。可視バンドから中間赤外バンドのように近赤外バンドを含む光学センサによって取得された衛星画像データには、大気と地形との双方の影響に基づく変動が含まれている。特に、起伏の大きい地域の可視バンド衛星画像データには大気状態の空間的変動に起因するパスラディアンスが大きく影響しているため、この影響を除去することは重要である。
【０００４】
非特許文献１では、近赤外バンドにおける各画素は可視バンドにおけるある画素に対応するものと仮定し、散乱が可視バンドで生じている場合は近赤外バンドにおける各画素は可視バンドにおける複数の画素に対応しているものと考えている。そこで、近赤外バンドで土地被覆分類等による分類により同じ分類クラスに分類された可視バンドの各画素の画素値を、それらの画素値の平均値で置き換えるという方法を大気補正法として提案している。この大気補正法は局所的にヘイズがある衛星画像データに対しては見かけ上良くなるが、局所的な情報は失われることになるという問題があった。例えば、近赤外バンドではヘイズのある地域での水田もクリアな地域での水田も同じ分類クラスに属しているが、可視バンドでは両者は異なる分類クラスに属している。ここで水田がヘイズのある地域に偏在していると、可視バンドにおける画素値の平均値はヘイズのある地域側に偏った値となる。このため、クリアな地域の水田の画素値を当該平均値で置き換えると、クリアな地域の水田も汚されてしまうことになる。つまり、非特許文献１の大気補正法では水田がヘイズのある地域に偏在しているというような局所的な情報は失われてしまうことになるという問題があった。さらに、非特許文献１では起伏のある地域を扱っていないため、起伏による陰影の情報という局所的な情報も失われてしまい、見かけも良くなくなってしまうという問題があった。
【０００５】
非特許文献２では、起伏の激しい山岳地域では地表面の傾斜による太陽入射照度の違いが重要な問題となる点に着目し、大気・地形効果補正の簡便な方法を提案している。詳しくは、パスラディアンスが標高の一次関数として表すことができるということを放射伝達のシミュレーションと実際の衛星画像データとを用いて示している。しかし、大気状態が水平方向に変化する場合は取り扱っていないという問題があった。
【０００６】
非特許文献３では、タセルドキャップ（Tasseled Cap）というランドサットＴＭの上述の各バンドの画素値を線形変換した指標の中の４番目の指標（簡便にはバンド１とバンド３との差）がヘイズの量と相関が高いという経験式を用いた方法が記載されている。しかし、非特許文献３では画素毎に分光反射特性が異なるということが無視されている。統計的にも変換後のバンド間の相関が高くなることが指摘されており、非特許文献３の方法は物理的に信頼性が低いという問題があった。
【０００７】
非特許文献４では、まず赤外バンドを分類し、その分類クラス毎にエアロゾル（ヘイズ）の光学的厚さを求め、次に、当該光学的厚さを空間的に平均した値を求めている。詳しくは、光学的厚さの空間的な平滑化を比較的狭い範囲（１５０ｍ程度）で通常の平均値を用いて行なっている。この平均値を利用して、画素値から地表面の反射率を補正している。しかし、非特許文献４では起伏のある地域を取り扱っていないという問題があった。さらに、光学的厚さおよび反射率という物理量を用いているため、クリアな地域における大気の散乱等を予め計算しておく必要がある。従って、クリアな地域を厳密に定めておく必要があり、煩雑で非実用的であるという問題があった。加えて、補正に利用する値は比較的狭い範囲における通常の平均値であるため、大きな範囲で平滑化を行った場合、必要とする変化を捉えることができないという問題があった。
【０００８】
【非特許文献１】M.J.Carlooto:“Reducingthe effects of space-varying, wavelength-dependent scattering in multispectralimagery”, Int.J.Remote Sensing, 1999, Vol.20, No.17, pp.3333-3344.
【非特許文献２】飯倉善和、横山隆三、「ランドサットＴＭ画像の大気および地形効果の補正」、日本リモートセンシング学会誌、１９９９年、第１９巻、第１号、ｐ．２−１６．(Y.Iikura and R.Yokoyama:“Correction of Atmospheric and TopographicEffects on Landsat TM Imagery”, Journal of Remote Sensing Society of Japan, 1999,Vol.19, No.1, pp.2-16.)
【非特許文献３】J.Lavreau:“De-hazingLandsat Thematic Mapper images”, Photogrammetric Engineering & RemoteSensing, 1991, Vol.57, No.10, pp.1297-1302.
【非特許文献４】S.Liang,H.Fang, and M.Chen:“Atmospheric Correction of Landsat ETM+ Land Surface Imagery- Part I：Methods”, IEEE Trans. on Geoscience and Remote Sensing, 2001, Vol.39,No.11, pp.2490-2498.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
上述のように、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データには大気状態の空間的変動に起因するパスラディアンスが大きく影響しているため、この影響を除去することは重要である。しかし、非特許文献１に記載された大気補正法では、局所的にヘイズがある衛星画像データに対しては見かけ上良くなるが、局所的な情報は失われることになるという問題があった。さらに、非特許文献１では起伏のある地域を扱っていないため、起伏による陰影の情報という局所的な情報も失われてしまい、見かけも良くなくなってしまうという問題があった。非特許文献２に記載された大気・地形効果補正の方法では、起伏のある地域を考慮しているが、大気状態が水平方向に変化する場合は取り扱っていないという問題があった。非特許文献３の方法は物理的に信頼性が低いという問題があった。非特許文献４に記載された方法では、起伏のある地域を取り扱っていないという問題があった。さらに、光学的厚さおよび反射率という物理量を用いているため、クリアな地域における大気の散乱等を予め計算しておく必要がある。従って、クリアな地域を厳密に定めておく必要があり、煩雑で非実用的であるという問題があった。加えて、補正に利用する値は比較的狭い範囲における通常の平均値であるため、大きな範囲で平滑化を行った場合、必要とする変化を捉えることができないという問題があった。
【００１０】
そこで、本発明の目的は、上記問題を解決するためになされたものであり、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することにある。さらに、本発明の目的は、簡便で実用的であり、大きな範囲で平滑化を行った場合でも、必要とする変化を捉えることができる画素値補正プログラム等を提供することにある。
【課題を解決するための手段】
【００１１】
この発明の画素値補正プログラムは、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムであって、コンピュータに、前記対象地域に関する近赤外バンドの衛星画像データを入力し、該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ、前記対象地域に関する可視バンドの衛星画像データを入力し、前記分類ステップで同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ、前記対象地域に関する数値標高モデルの座標データを入力し、前記分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ、前記分類クラス毎に、前記可視バンドデータ抽出ステップで抽出された画素値、前記数値標高データ抽出ステップで抽出された標高及び所定の基準標高に基づき画素毎に大気影響度を求める分析ステップ、前記分析ステップで求めた各画素の大気影響度に対して中央値フィルタを適用し、該画素毎に地域的な大気影響度を求めるフィルタステップ、前記可視バンドデータ抽出ステップで抽出した各画素の画素値と、前記数値標高データ抽出ステップで抽出した標高と、前記所定の基準標高と、前記フィルタステップで求めた該画素の地域的な大気影響度とに基づき、各画素の画素値を補正する補正ステップを実行させるための画素値補正プログラムである。
【００１２】
ここで、この発明の画素値補正プログラムにおいて、前記所定の基準は、近赤外バンドの衛星画像データの画素値が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準とすることができる。
【００１３】
ここで、この発明の画素値補正プログラムにおいて、前記所定の基準標高は前記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値とすることができる。
【００１４】
ここで、この発明の画素値補正プログラムにおいて、前記分析ステップは、前記分類クラス毎に、横軸を標高、縦軸を画素値とするグラフを用い各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットするグラフ作成ステップと、クリアな領域における傾きＢ０（画素値／標高）を設定する傾き設定ステップと、前記傾き設定ステップで設定された傾きＢ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算する仮の画素値計算ステップと、前記仮の画素値計算ステップで計算された仮の画素値Ｙｐｉの中央値を求め該中央値を基準標高Ｘ０における画素値Ｙ０とする画素値決定ステップと、前記画素値決定ステップで決定された画素値Ｙ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算し該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とするステップとを備えることができる。
【００１５】
ここで、この発明の画素値補正プログラムにおいて、前記フィルタステップは、各画素の周囲６００ｍ程度以内の画素の大気影響度を抽出し、該大気影響度の中央値を用いて平滑化した結果を各画素の地域的な大気影響度として求めることができる。
【００１６】
ここで、この発明の記載の画素値補正プログラムにおいて、前記補正ステップは、各画素の標高Ｘｉ、画素値Ｙｉ、基準標高Ｘ０及び地域的な大気影響度Ｂｉ’を用いて補正後の画素値Ｙｉ’を、Ｙｉ’＝Ｙｉ―Ｂｉ’×（Ｘ０−Ｘｉ）、ｉ＝１〜全画素数、により計算することができる。
【００１７】
この発明の記録媒体は、本発明の画素値補正プログラムを記録したコンピュータ読取り可能な記録媒体である。
【発明の効果】
【００１８】
本発明の画素値補正プログラムによれば、大気影響度を土地被覆分類等の分類クラス毎に標高を加味した式により計算しているため、タセルドキャップよりも物理的に信頼性の高い画素値補正プログラム等を提供することができる。画素毎の大気影響度から地域的な大気影響度を求める際に中央値フィルタを利用しているため、ヘイズ等の非連続な大気の状況の変化を自然に把握することができ、大きい範囲で平滑化を行った場合であっても、必要な変化を捉えることができる。さらに、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めた地域的な大気影響度を求めることができる。画素における大気影響度または地域的な大気影響度を基準標高における画素に対する傾きと考えたことにより、補正量を標高の一次式として簡単に求めることができる。
【００１９】
この結果、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することができるという効果がある。さらに、光学的厚さまたは反射率等の物理量を用いていないため、簡便で実用的な画素値補正プログラム等を提供することができるという効果がある。
【発明を実施するための最良の形態】
【００２０】
まず、本発明の、対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムの基盤となる理論について説明し、次に各実施例について図面を参照して詳細に説明する。
【００２１】
理論の説明．
１．背景技術の非特許文献２で説明した、パスラディアンスを標高の一次関数として表すことができるという点を拡張して、当該一次関数の傾きが場所により異なるという以下の式1および２を本発明の画素値補正プログラムの基本モデルとした。これにより、大気の補正を見通しのよいものとした。標高（ｍ）をＸ軸、画素値（Digital Number :ＤＮ）をＹ軸とし、標高Ｘｉにおける補正前の画素値をＹｉ、補正量をＤｉ、傾きをＢｉ、パラメータをＡ、補正された画素値をＹｉ’とすると、
【００２２】
Ｙｉ’＝Ｙｉ−Ｄｉ （１）
【００２３】
と表すことができる。補正量Ｄｉは、
【００２４】
Ｄｉ＝Ａ−Ｂｉ×Ｘｉ （２）
【００２５】
である。ここで、大気の影響が無視できる標高（基準標高Ｘ０）での補正量を０とすれば、Ａは以下のように設定できる。
【００２６】
０＝Ａ−Ｂｉ×Ｘ０
∴Ａ＝Ｂｉ×Ｘ０
【００２７】
上記パラメータＡを式１に代入すると、式１は、以下の式３のようになり、パラメータはＢｉのみでよいことになる。
【００２８】
Ｄｉ＝Ｂｉ×Ｘ０−Ｂｉ×Ｘｉ＝Ｂｉ×（Ｘ０−Ｘｉ） （３）
【００２９】
対象画素（Ｘｉ，Ｙｉ）の傾きＢｉは、基準標高Ｘ０における画素値Ｙ０（大気の影響がない場合の画素値）がわかれば、以下の式４から求めることができる。
【００３０】
Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ） （４）
【００３１】
２．基本モデルのＢｉを画素（Ｘｉ，Ｙｉ）毎に求めるためには、式４の計算を行うことになる。ここで、式４における画素値Ｙ０は土地被覆分類等により分類された分類クラス毎に必要である。本発明の画素値補正プログラムでは、画素値Ｙ０を求めるために、背景技術の非特許文献１で説明したように近赤外バンドで同じ分類クラスに分類された可視バンドの各画素の画素値を利用する。非特許文献１では可視バンドの各画素の平均値に各画素値を置き換えているが、平均値からの差が大気の影響度と考えることもできる。画像の半分以上がクリアな領域である場合、平均値ではなくそれらの中央値に画素値を置き換えた方が良い。画素値Ｙ０を求めるためにクリアな領域での傾きＢ０を利用する。Ｂ０は６Ｓ（Second Simulation of the Satellite Signal in the Solar Spectrumの略で、フランスのLille科学技術大学で開発されたプログラム。）等により設定した値（バンド１で７〜８ＤＮ／ｋｍ）を用いる。まず、同じ分類クラスに含まれるすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを以下の式５を用いて計算する。
【００３２】
Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ） （５）
【００３３】
クリアな領域に含まれる画素のＹｐｉは標高によらず一定となるため、これを大気の影響がない場合における画素値Ｙ０とみなすことができる。分類クラスの中でクリアな領域に含まれる画素が多い場合、Ｙｐｉの中央値はＹ０となる。クリアでない領域に含まれる画素のＹｐｉは大きくなるため、このような画素が多い場合には、すべての画素（Ｘｉ，Ｙｉ）について計算されたＹｐｉの中央値よりもクリアな地域について計算された中央値を使用すればよい。具体的には、クリアな地域が占める画像中の割合をＷ（％）とした場合、すべての画素（Ｘｉ，Ｙｉ）について計算されたＹｐｉの度数分布の下からＷ／２（％）のＹｐｉを用いて計算された中央値を使用すればよい。
【００３４】
３．上述した１および２により画素毎に傾きＢｉを求めることができる。傾きＢｉは大気の影響度として採用することができる。但し、傾きＢｉには種々の誤差要因が含まれている。衛星画像データからは、大気の影響度には地域的な特徴があり、ある場所ではクリアであるが他の場所では薄い霧またはヘイズがかかっていることがわかる。大気の影響度はある程度の広がりで緩やかに変化すると考えられる。従って、大気の影響度としては画素単位ではなくその近傍のデータも利用した方が誤差を少なくすることができると考えられる。そこで、対象画素の周囲の画素の傾きＢｉを抽出してそれらの中央値を利用するという空間的な中央値フィルタを適用して、地域的な大気影響度Ｂｉ’を求める。平均値ではなく中央値を利用する理由は、クリアな領域とクリアでない領域とがはっきりと分かれている場合、中央値に比べて平均値を利用すると両領域の境界がぼけてしまうためである。この地域的な大気影響度Ｂｉ’を用いて以下の式６により補正量Ｄｉを求める。
【００３５】
Ｄｉ＝Ｂｉ’×（Ｘ０−Ｘｉ） （６）
【００３６】
補正された画素値Ｙｉ’は、式１および６により、以下の式７のようになる。
【００３７】
Ｙｉ’＝Ｙｉ−Ｂｉ’×（Ｘ０−Ｘｉ） （７）
【実施例１】
【００３８】
図１は、本発明の実施例１における画素値補正プログラムの流れをフローチャートで示す。以下では図１における各処理ブロックの説明を、図２〜４の実例、図５に示される細部のフローチャート、図６に示されるグラフ、図７〜９の実例を用いて説明する。
【００３９】
図１に示されるように、コンピュータに、まず、所望の対象地域（例えば岩手県大川地区）に関する近赤外バンドの衛星画像データを入力し、当該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ（ステップＳ１０）を実行させる。入力する近赤外バンドの衛星画像データは一般的なものであるため、説明は省略する。所定の基準としては、例えば近赤外バンドの衛星画像データの画素値（ＤＮ）が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準を用いることができる。他の基準を用いてもよいことは勿論である。
【００４０】
図２は、分類ステップ（ステップＳ１０）を実行させた結果の分類データ画像を示す。図２に示される分類データ画像はカラー合成画像である。
【００４１】
次に、図１に示されるように、コンピュータに、上記対象地域と同じ対象地域に関する可視バンドの衛星画像データを入力し、分類ステップ（ステップＳ１０）で同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ（ステップＳ２０）を実行させる。
【００４２】
図３は、可視バンド分類ステップ（ステップＳ２０）で入力する可視バンド（バンド１）の衛星画像データの一例を示す。図３の右上部Ｈａｚｅで示されるように、当該部分はヘイズの影響により明るくなっている。
【００４３】
続いて、図１に示されるように、コンピュータに、上記対象地域と同じ対象地域に関する数値標高モデルの座標データを入力し、分類ステップ（ステップＳ１０）で分類された分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ（ステップＳ３０）を実行させる。
【００４４】
図４は、画素値取得ステップ（ステップＳ３０）で入力する数値標高モデルの座標データ画像の一例を示す。当該座標データ画像は、国土地理院の数値データ５０ｍメッシュ（標高）からＵＴＭ座標系（Universal Transverse Mercator system）へ変換して作成したものである。
【００４５】
次に、図１に示されるように、コンピュータに、上記分類クラス毎に、可視バンドデータ抽出ステップ（ステップＳ２０）で抽出された画素値と、数値標高データ抽出ステップ（ステップＳ３０）で抽出された標高と、所定の基準標高とに基づき画素毎に大気影響度を求める分析ステップ（ステップＳ４０）を実行させる。
【００４６】
図５は分析ステップ（ステップＳ４０）の詳細な流れをフローチャートで示し、図６は分析ステップ（ステップＳ４０）の詳細を説明するためのグラフを示す。以下の処理は上記分類クラス毎に実行される。図６に示される横軸を標高（ｍ）、縦軸を画素値（ＤＮ）とするグラフを用い、各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットする（グラフ作成ステップ。ステップＳ４１）。
【００４７】
次に、クリアな領域における傾きＢ０（画素値／標高）を設定する（傾き設定ステップ。ステップＳ４２）。上述のように、Ｂ０は６Ｓ等の放射伝達シミュレーションにより設定した値（バンド１で７〜８ＤＮ／ｋｍ）を用いればよい。
【００４８】
傾き設定ステップ（ステップＳ４２）で設定された傾きＢ０を用い、グラフ作成ステップ（ステップＳ４１）でプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、上述の式５、
【００４９】
Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜上記分類クラスの画素数 （５）
【００５０】
を用いて計算する（仮の画素値計算ステップ。ステップＳ４３）。所定の基準標高Ｘ０としては、上記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値とすることが好適である。図６では所定の基準標高Ｘ０として１２００ｍを用いているが、所定の基準標高Ｘ０は１２００ｍに限定されるものではない。
【００５１】
仮の画素値計算ステップ（ステップＳ４３）で計算された仮の画素値Ｙｐｉの中央値を求め、当該中央値を基準標高Ｘ０における画素値Ｙ０とする（画素値決定ステップ。ステップＳ４４）。図６では、画素値Ｙ０は約７５（ＤＮ）となっている。
【００５２】
画素値決定ステップ（ステップＳ４４）で決定された画素値Ｙ０を用い、グラフ作成ステップ（ステップＳ４１）でプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、上述の式４、
【００５３】
Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜上記分類クラスの画素数 （４）
【００５４】
を用いて計算し、当該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とする（ステップＳ４５）。図６では画素（４５０、９０）について基準標高１２００ｍにおける画素（１２００、７５）に対する傾きＢｉが直線ｉにより示されている。この傾きＢｉ、即ち７５０ｍ（＝１２００ｍ−４５０ｍ）当たり１５ＤＮ（＝９０ＤＮ−７５ＤＮ）変わるという値が画素（４５０、９０）における大気影響度である。図７は、上述のようにして求められた画素毎の大気影響度Ｂｉを示す画像である。
【００５５】
図１に戻り、続いてコンピュータに、分析ステップ（ステップＳ４０）で求めた各画素（Ｘｉ，Ｙｉ）の大気影響度Ｂｉに対して中央値フィルタを適用し、画素（Ｘｉ，Ｙｉ）毎に地域的な大気影響度Ｂｉ’を求めるフィルタステップ（ステップＳ５０）を実行させる。図８は、フィルタステップ（ステップＳ５０）の実行結果を示す画像である。図８では、フィルタステップ（ステップＳ５０）において、各画素（Ｘｉ，Ｙｉ）の周囲６００ｍ程度以内の画素の大気影響度Ｂｉを抽出し、この大気影響度Ｂｉの中央値を用いて平滑化した結果を各画素（Ｘｉ，Ｙｉ）の地域的な大気影響度Ｂｉ’として求めている。
【００５６】
最後に、コンピュータに、可視バンドデータ抽出ステップ（ステップ２０）で抽出した各画素の画素値Ｙｉと、数値標高データ抽出ステップ（ステップＳ３０）で抽出した標高Ｘｉと、所定の基準標高Ｘ０と、フィルタステップ（ステップ５０）で求めた当該画素の地域的な大気影響度Ｂｉ’とに基づき、上述の式６、
【００５７】
Ｄｉ＝Ｂｉ’×（Ｘ０−Ｘｉ） （６）
【００５８】
により補正量Ｄｉを求め、上述の式７、
【００５９】
Ｙｉ’＝Ｙｉ−Ｂｉ’×（Ｘ０−Ｘｉ） （７）
【００６０】
により各画素の画素値Ｙｉを補正してＹｉ’を求める補正ステップ（ステップ６０）を実行させる。ここで、ｉ＝１〜全画素数である。図９は、上述した補正ステップ（ステップＳ６０）を実行した結果を示す画像である。
【００６１】
以上より、本発明の実施例１によれば、コンピュータに、まず、所望の対象地域に関する近赤外バンドの衛星画像データを入力し、当該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップを実行させる。次に、上記対象地域と同じ対象地域に関する可視バンドの衛星画像データを入力し、分類ステップ（ステップＳ１０）で同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ（ステップＳ２０）を実行させる。続いて、上記対象地域と同じ対象地域に関する数値標高モデルの座標データを入力し、上記分類クラス毎に、各画素の標高Ｘｉを抽出する数値標高データ抽出ステップ（ステップＳ３０）を実行させる。上記分類クラス毎に、可視バンドデータ抽出ステップ（ステップＳ３０）で抽出された画素値と、数値標高データ抽出ステップ（ステップＳ３０）で抽出された標高と、所定の基準標高とに基づき画素毎に大気影響度を求める分析ステップ（ステップＳ４０）を実行させる。分析ステップ（ステップＳ４０）で求めた各画素（Ｘｉ，Ｙｉ）の大気影響度Ｂｉに対して中央値フィルタを適用し、画素（Ｘｉ，Ｙｉ）毎に地域的な大気影響度Ｂｉ’を求めるフィルタステップ（ステップＳ５０）を実行させる。可視バンドデータ抽出ステップ（ステップＳ２０）で抽出した各画素の画素値Ｙｉと、数値標高データ抽出ステップ（ステップＳ３０）で抽出した標高と、所定の基準標高Ｘ０と、フィルタステップ（ステップ５０）で求めた当該画素の地域的な大気影響度Ｂｉ’とに基づき、上述の式６により補正量Ｄｉを求め、上述の式７により各がその画素値Ｙｉを補正してＹｉ’を求める補正ステップ（ステップＳ６０）を実行させる。
【００６２】
本発明の画素値補正プログラムでは、大気影響度Ｂｉを土地被覆分類等の分類クラス毎に標高Ｘｉを加味した式４により計算しているため、タセルドキャップよりも物理的に信頼性の高い画素値補正プログラム等を提供することができる。本発明の画素値補正プログラムでは、画素毎の大気影響度Ｂｉから地域的な大気影響度Ｂｉ’を求める際に中央値フィルタを利用しているため、ヘイズ等の非連続な大気の状況の変化を自然に把握することができ、大きい範囲で平滑化を行った場合であっても、必要な変化を捉えることができる。さらに、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めた地域的な大気影響度Ｂｉ’を求めることができる。本発明の画素値補正プログラムでは、画素（Ｘｉ，Ｙｉ）における大気影響度Ｂｉまたは地域的な大気影響度Ｂｉ’を基準標高Ｘ０における画素（Ｘ０，Ｙ０）に対する傾きＢｉと考えたことにより、補正量Ｄｉを標高Ｘｉの一次式として簡単に求めることができる。以上より、近赤外バンドを含む光学センサによって取得された衛星画像データ、特に起伏の大きい地域の可視バンド衛星画像データから、局所的な情報を失なわず見かけが悪くなることがなく大気状態が水平方向に変化する場合も含めてパスラディアンスの影響を除去することができる物理的に信頼性の高い画素値補正プログラム等を提供することができる。さらに、光学的厚さまたは反射率等の物理量を用いていないため、簡便で実用的な画素値補正プログラム等を提供することができる。
【実施例２】
【００６３】
図１０は、上述した実施例１を実現するための本発明のコンピュータ・プログラム（画素値補正プログラム）を実行するコンピュータの内部回路１０を示すブロック図である。図１０に示されるように、ＣＰＵ１１、ＲＯＭ１２、ＲＡＭ１３、画像制御部１４、コントローラ１８、入力制御部２０および外部インタフェース（Ｉ／Ｆ）部２１はバス２２に接続されている。図１０において、上述の本発明の画素値補正プログラムは、ＲＯＭ１２、ＨＤＤ１７ａまたはＤＶＤ１７ｎ（あるいはＣＤ−ＲＯＭ）等の記録媒体（脱着可能な記録媒体を含む）に記録されている。ＨＤＤ１７ａまたはＤＶＤ１７ｎ等には入力した近赤外バンドの衛星画像データ、可視バンドの衛星画像データ、数値標高モデルの座標データ等が記録されている。本発明の画素値補正プログラムは、ＲＯＭ１２からバス２２を介し、あるいはＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体からコントローラ１８を経由してバス２２を介しＲＡＭ１３へロードされる。入力制御部２０はマウス、テンキー等の入力操作部１９と接続され入力制御等を行う。画像メモリであるＶＲＡＭ１５は、分類データ画像、画素毎の大気影響度Ｂｉを示す画像、フィルタステップ（ステップＳ５０）の実行結果を示す画像、補正ステップ（ステップＳ６０）を実行した結果を示す画像等を表示するディスプレイ等の表示部１６の少なくとも一画面分のデータ容量に相当する容量を有しており、画像制御部１６はＶＲＡＭ１５のデータを画像データへ変換して表示部１６へ送出する機能を有している。外部Ｉ／Ｆ部２１は、近赤外バンドの衛星画像データ、可視バンドの衛星画像データ、数値標高モデルの座標データ等を入力する際におけるインタフェース機能を有する。
【００６４】
上述のようにＣＰＵ１１が上述の本発明の画素値補正プログラムを実行することにより、本発明の目的を達成することができる。当該画素値補正プログラムは上述のようにＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体の形態でコンピュータＣＰＵ１１に供給することができ、当該画素値補正プログラムを記録したＨＤＤ１７ａまたはＤＶＤ１７ｎ等の記録媒体も同様に本発明を構成することになる。当該画素値補正プログラムを記録した記録媒体としては上述された記録媒体の他に、例えばメモリ・カード、メモリ・スティック等を用いることができる。
【産業上の利用可能性】
【００６５】
本発明の活用例として、ランドサットＴＭセンサにより取得された衛星画像データの画素値の補正に適用することができる。
【図面の簡単な説明】
【００６６】
【図１】本発明の実施例１における補正プログラムの流れを示すフローチャートである。
【図２】分類ステップ（ステップＳ１０）を実行させた結果の分類データ画像を示す図である。
【図３】可視バンド分類ステップ（ステップＳ２０）で入力する可視バンド（バンド１）の衛星画像データの一例を示す図である。
【図４】分析ステップ（ステップＳ４０）で入力する数値標高モデルの座標データ画像の一例を示す図である。
【図５】分析ステップ（ステップＳ４０）の詳細な流れを示すフローチャートである。
【図６】分析ステップ（ステップＳ４０）の詳細を説明するためのグラフである。
【図７】分析ステップ（ステップＳ４０）で求められた画素毎の大気影響度Ｂｉを示す画像である。
【図８】フィルタステップ（ステップＳ５０）の実行結果を示す画像である。
【図９】補正ステップ（ステップＳ６０）を実行した結果を示す画像である。
【図１０】本発明の実施例を実現するための本発明の画素値補正プログラムを実行するコンピュータの内部回路１０を示すブロック図である。
【符号の説明】
【００６７】
Ｂｉ、Ｂ０ 傾き、 Ｈａｚｅ ヘイズによる影響がある領域、 Ｘ０ 基準標高、 Ｙ０ 基準標高における画素値、 １０ 内部回路、 １１ ＣＰＵ、 １２ ＲＯＭ、 １３ ＲＡＭ、 １４ 画像制御部、 １５ ＶＲＡＭ、 １６ 表示部、 １７ａ ＨＤＤ、 １７ｎ ＤＶＤ、 １８ コントローラ、 １９ 入力操作部、 ２０ 入力制御部、 ２１ 外部Ｉ／Ｆ部、 ２２ バス。
【特許請求の範囲】
【請求項１】
対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムであって、コンピュータに、
前記対象地域に関する近赤外バンドの衛星画像データを入力し、該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ、
前記対象地域に関する可視バンドの衛星画像データを入力し、前記分類ステップで同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ、
前記対象地域に関する数値標高モデルの座標データを入力し、前記分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ、
前記分類クラス毎に、前記可視バンドデータ抽出ステップで抽出された画素値、前記数値標高データ抽出ステップで抽出された標高及び所定の基準標高に基づき画素毎に大気影響度を求める分析ステップ、
前記分析ステップで求めた各画素の大気影響度に対して中央値フィルタを適用し、該画素毎に地域的な大気影響度を求めるフィルタステップ、
前記可視バンドデータ抽出ステップで抽出した各画素の画素値と、前記数値標高データ抽出ステップで抽出した標高と、前記所定の基準標高と、前記フィルタステップで求めた該画素の地域的な大気影響度とに基づき、各画素の画素値を補正する補正ステップを実行させるための画素値補正プログラム。
【請求項２】
請求項１記載の画素値補正プログラムにおいて、前記所定の基準は、近赤外バンドの衛星画像データの画素値が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準であることを特徴とする画素値補正プログラム。
【請求項３】
請求項１又は２記載の画素値補正プログラムにおいて、前記所定の基準標高は前記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値であることを特徴とする画素値補正プログラム。
【請求項４】
請求項１ないし３のいずれかに記載の画素値補正プログラムにおいて、前記分析ステップは、前記分類クラス毎に、
横軸を標高、縦軸を画素値とするグラフを用い各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットするグラフ作成ステップと、
クリアな領域における傾きＢ０（画素値／標高）を設定する傾き設定ステップと、
前記傾き設定ステップで設定された傾きＢ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算する仮の画素値計算ステップと、
前記仮の画素値計算ステップで計算された仮の画素値Ｙｐｉの中央値を求め該中央値を基準標高Ｘ０における画素値Ｙ０とする画素値決定ステップと、
前記画素値決定ステップで決定された画素値Ｙ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算し該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とするステップとを備えたことを特徴とする画素値補正プログラム。
【請求項５】
請求項１ないし４のいずれかに記載の画素値補正プログラムにおいて、前記フィルタステップは、各画素の周囲６００ｍ程度以内の画素の大気影響度を抽出し、該大気影響度の中央値を用いて平滑化した結果を各画素の地域的な大気影響度として求めることを特徴とする画素値補正プログラム。
【請求項６】
請求項１ないし５のいずれかに記載の画素値補正プログラムにおいて、前記補正ステップは、各画素の標高Ｘｉ、画素値Ｙｉ、基準標高Ｘ０及び地域的な大気影響度Ｂｉ’を用いて補正後の画素値Ｙｉ’を、Ｙｉ’＝Ｙｉ―Ｂｉ’×（Ｘ０−Ｘｉ）、ｉ＝１〜全画素数、により計算することを特徴とする画素値補正プログラム。
【請求項７】
請求項１ないし６のいずれかに記載の画素値補正プログラムを記録したコンピュータ読取り可能な記録媒体。
【請求項１】
対象地域に関する可視バンドの衛星画像データの画素値を補正するための画素値補正プログラムであって、コンピュータに、
前記対象地域に関する近赤外バンドの衛星画像データを入力し、該近赤外バンドの衛星画像データの各画素を所定の基準に基づき分類クラスに分類する分類ステップ、
前記対象地域に関する可視バンドの衛星画像データを入力し、前記分類ステップで同じ分類クラスに分類された各画素に対応する可視バンドの衛星画像データの画素値を抽出する可視バンドデータ抽出ステップ、
前記対象地域に関する数値標高モデルの座標データを入力し、前記分類クラス毎に各画素の標高を抽出する数値標高データ抽出ステップ、
前記分類クラス毎に、前記可視バンドデータ抽出ステップで抽出された画素値、前記数値標高データ抽出ステップで抽出された標高及び所定の基準標高に基づき画素毎に大気影響度を求める分析ステップ、
前記分析ステップで求めた各画素の大気影響度に対して中央値フィルタを適用し、該画素毎に地域的な大気影響度を求めるフィルタステップ、
前記可視バンドデータ抽出ステップで抽出した各画素の画素値と、前記数値標高データ抽出ステップで抽出した標高と、前記所定の基準標高と、前記フィルタステップで求めた該画素の地域的な大気影響度とに基づき、各画素の画素値を補正する補正ステップを実行させるための画素値補正プログラム。
【請求項２】
請求項１記載の画素値補正プログラムにおいて、前記所定の基準は、近赤外バンドの衛星画像データの画素値が一致しているか又は所定の範囲内にある画素を同じ分類クラスに分類する基準であることを特徴とする画素値補正プログラム。
【請求項３】
請求項１又は２記載の画素値補正プログラムにおいて、前記所定の基準標高は前記分類クラスの画素の標高より高く１０００ｍから１５００ｍ程度の範囲内の値であることを特徴とする画素値補正プログラム。
【請求項４】
請求項１ないし３のいずれかに記載の画素値補正プログラムにおいて、前記分析ステップは、前記分類クラス毎に、
横軸を標高、縦軸を画素値とするグラフを用い各画素の標高Ｘｉ及び画素値Ｙｉに対応する座標（Ｘｉ，Ｙｉ）に画素をプロットするグラフ作成ステップと、
クリアな領域における傾きＢ０（画素値／標高）を設定する傾き設定ステップと、
前記傾き設定ステップで設定された傾きＢ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高Ｘ０における仮の画素値Ｙｐｉを、Ｙｐｉ＝Ｙｉ−Ｂ０×（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算する仮の画素値計算ステップと、
前記仮の画素値計算ステップで計算された仮の画素値Ｙｐｉの中央値を求め該中央値を基準標高Ｘ０における画素値Ｙ０とする画素値決定ステップと、
前記画素値決定ステップで決定された画素値Ｙ０を用い、前記グラフ作成ステップでプロットされたすべての画素（Ｘｉ，Ｙｉ）について、基準標高における画素（Ｘ０，Ｙ０）に対する傾きＢｉを、Ｂｉ＝（Ｙｉ−Ｙ０）／（Ｘ０−Ｘｉ）、ｉ＝１〜前記分類クラスの画素数、を用いて計算し該傾きＢｉを画素（Ｘｉ，Ｙｉ）の大気影響度とするステップとを備えたことを特徴とする画素値補正プログラム。
【請求項５】
請求項１ないし４のいずれかに記載の画素値補正プログラムにおいて、前記フィルタステップは、各画素の周囲６００ｍ程度以内の画素の大気影響度を抽出し、該大気影響度の中央値を用いて平滑化した結果を各画素の地域的な大気影響度として求めることを特徴とする画素値補正プログラム。
【請求項６】
請求項１ないし５のいずれかに記載の画素値補正プログラムにおいて、前記補正ステップは、各画素の標高Ｘｉ、画素値Ｙｉ、基準標高Ｘ０及び地域的な大気影響度Ｂｉ’を用いて補正後の画素値Ｙｉ’を、Ｙｉ’＝Ｙｉ―Ｂｉ’×（Ｘ０−Ｘｉ）、ｉ＝１〜全画素数、により計算することを特徴とする画素値補正プログラム。
【請求項７】
請求項１ないし６のいずれかに記載の画素値補正プログラムを記録したコンピュータ読取り可能な記録媒体。
【図１】
【図５】
【図１０】
【図２】
【図３】
【図４】
【図６】
【図７】
【図８】
【図９】
【図５】
【図１０】
【図２】
【図３】
【図４】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２００７−１７２５９０（Ｐ２００７−１７２５９０Ａ）
【公開日】平成１９年７月５日（２００７．７．５）
【国際特許分類】
【出願番号】特願２００６−３０９００９（Ｐ２００６−３０９００９）
【出願日】平成１８年１１月１５日（２００６．１１．１５）
【出願人】（５０４２２９２８４）国立大学法人弘前大学 (162)
【出願人】（５９１１４６２３９）いであ株式会社 (11)
【Ｆターム（参考）】
【公開日】平成１９年７月５日（２００７．７．５）
【国際特許分類】
【出願日】平成１８年１１月１５日（２００６．１１．１５）
【出願人】（５０４２２９２８４）国立大学法人弘前大学 (162)
【出願人】（５９１１４６２３９）いであ株式会社 (11)
【Ｆターム（参考）】
[ Back to top ]