JP6463244B2 - 画像処理装置及び画像処理方法 - Google Patents

Description

本発明は、複数の観測波長帯の電磁波を観測して得られるマルチスペクトル画像を処理する画像処理技術に関するものである。

マルチスペクトルセンサやハイパースペクトルセンサと呼ばれる光学センサは、観測対象物から放射された電磁波のうち複数の観測波長帯（バンド）の電磁波を検出することにより、観測波長帯ごとの画像情報を有するマルチスペクトル画像を生成することができる。この種の光学センサは、人工衛星、航空機並びに宇宙探査機などに搭載され、地球などの惑星を観測する際に使用されている。また、この種の光学センサを用いて地球が観測される場合、地表における観測対象物の分光反射特性の違いに基づいて、当該観測対象物もしくは現象の識別、並びに当該観測対象物の種々の分類を行うことが可能である。地表から放射された電磁波は、大気を伝搬する過程で散乱や減衰などの影響を受けた後に光学センサにより観測されるので、その影響を低減させるためにマルチスペクトル画像を補正する必要がある。この種の画像補正技術は、たとえば、特許文献１（特開２００５−１５７５６１号公報）及び特許文献２（特開２０１３−２２５２４３号公報）に開示されている。

特許文献１の従来技術は、観測対象物である被写体の陰影部の放射輝度から大気伝搬補正パラメータを算出している。しかしながら、この従来技術は、ユーザが目視で陰影部を指定して、その陰影部の放射輝度を抽出することを前提とする技術である。このため、ユーザが被写体の陰影部を適正に指定しなければ、適正な大気伝搬補正パラメータを算出することができず、マルチスペクトル画像を適正に補正することができない。

一方、特許文献２の従来技術は、センサパラメータ記憶部に記憶されている既知の大気散乱放射輝度比α_ｉを用いて大気散乱放射輝度を算出し、この大気散乱放射輝度と既知の大気透過率とを用いて各画素の放射輝度値を自動的に補正することができる。ここで、大気散乱放射輝度は、大気中の分子または粒子によって散乱された太陽光がセンサに入射することで発生する放射輝度をいう。

特開２００５−１５７５６１号公報（たとえば、段落０００７，００１３） 特開２０１３−２２５２４３号公報（たとえば、段落００２１〜００２４）

しかしながら、特許文献２の従来技術では、大気散乱放射輝度比α_ｉは、観測波長帯に依存する成分に過ぎない。電磁波の大気伝搬に起因する分光反射特性の変化量は、エアロゾルと呼ばれる浮遊微粒子や気相分子を含む大気の状態に応じて異なるところ、特許文献２の従来技術では、霧や靄（もや）などによる太陽光の大気散乱の影響が観測画像内に一様に現れない場合には、観測画像を適正に補正することが難しいという課題がある。

上記に鑑みて本発明の目的は、霧や靄などによる太陽光の大気散乱の影響が観測画像内に空間的に一様に現れない場合でも、当該観測画像を適正に補正することができる画像処理装置を提供する点にある。

本発明の一態様による画像処理装置は、大気を有する惑星の地表の被撮像領域を観測して複数の観測波長帯にそれぞれ対応する複数のバンド画像を生成する光学センサから、当該複数のバンド画像の供給を受ける画像入力部と、前記複数のバンド画像のうちの少なくとも１つのバンド画像の画素値に基づいて、前記被撮像領域における暗部領域を検出する暗部検出部と、前記観測波長帯に依存する散乱特性データが記憶されているデータ記憶部と、前記複数のバンド画像における当該暗部領域の画素値に基づき、前記散乱特性データを用いて当該暗部領域に対応する大気の散乱度を推定する散乱度推定部と、当該推定された散乱度に基づいて、前記被撮像領域における当該暗部領域以外の他の領域に対応する大気の散乱度を補間することにより、前記被撮像領域全体に対応する大気の散乱度分布を生成する散乱度補間部と、前記散乱度分布を用いて前記複数のバンド画像の画素値を補正することにより、前記複数の観測波長帯にそれぞれ対応する複数の補正画像を算出する補正画像算出部とを備えることを特徴とする。

本発明の他の態様による画像処理方法は、大気を有する惑星の地表の被撮像領域を観測して複数の観測波長帯にそれぞれ対応する複数のバンド画像を生成する光学センサから当該複数のバンド画像の供給を受ける画像処理装置において実行される画像処理方法であって、前記複数のバンド画像のうちの少なくとも１つのバンド画像の画素値に基づいて前記被撮像領域における暗部領域を検出するステップと、前記惑星の大気散乱モデルに基づいて計算された、前記観測波長帯に依存する散乱特性データが記憶されているデータ記憶部から、前記散乱特性データを取得するステップと、前記複数のバンド画像における当該暗部領域の画素値に基づき、前記散乱特性データを用いて当該暗部領域に対応する大気の散乱度を推定するステップと、当該推定された散乱度に基づいて、前記被撮像領域における当該暗部領域以外の他の領域に対応する大気の散乱度を補間することにより、前記被撮像領域全体に対応する大気の散乱度分布を生成するステップと、前記散乱度分布を用いて前記複数のバンド画像の画素値を補正するステップとを備えることを特徴とする。

本発明によれば、暗部領域に対応する大気の散乱度が算出され、当該散乱度に基づいて被撮像領域全体に対応する散乱度分布が算出される。この散乱度分布を用いてバンド画像を補正することができる。したがって、太陽光の大気散乱の影響がバンド画像内に空間的に一様に現れない場合でも、当該バンド画像を適正に補正することができる。

本発明に係る実施の形態１の画像処理装置の概略構成を示す機能ブロック図である。 マルチスペクトル画像の構成を示す図である。 観測衛星によって地球の地表が観測される様子を模式的に示す図である。 赤色波長域、緑色波長域及び青色波長域のバンド画像のヒストグラムの例を示すグラフである。 正規分布とバンド画像のヒストグラムとの間の関係を示すグラフである。 ＮＩＲ（近赤外線）、赤色波長域、緑色波長域及び青色波長域のバンド画像のヒストグラムの例を示すグラフである。 大気散乱放射輝度の分布を示すグラフである。 実施の形態１に係る画像処理の手順の一例を概略的に示すフローチャートである。 実施の形態１の画像処理装置のハードウェア構成例を示す機能ブロック図である。 本発明に係る実施の形態２の画像処理装置の概略構成を示す機能ブロック図である。 実施の形態２に係る画像処理の手順の一例を概略的に示すフローチャートである。 本発明に係る実施の形態３の画像処理装置の概略構成を示す機能ブロック図である。 実施の形態３に係る画像処理の手順の一例を概略的に示すフローチャートである。

以下、図面を参照しつつ、本発明に係る種々の実施の形態について詳細に説明する。なお、図面において同一符号を付された構成要素は、同一機能及び同一構成を有するものとする。

実施の形態１．
図１は、本発明に係る実施の形態１の画像処理装置１の概略構成を示す機能ブロック図である。この画像処理装置１は、マルチスペクトルセンサまたはハイパースペクトルセンサと呼ばれる光学センサからマルチスペクトル画像の供給を受ける画像入力部１０を備えている。光学センサは、大気を有する惑星の地表の被撮像領域を観測することができるように人工衛星、航空機または宇宙探査機などのプラットフォームに搭載されている。画像入力部１０は、この光学センサを搭載するプラットフォームと通信を行う機能を有し、当該プラットフォームからマルチスペクトル画像ＭＳＩのデータを受信することができる。

図２は、マルチスペクトル画像ＭＳＩの構成例を概略的に示す図である。図２に示されるように、マルチスペクトル画像ＭＳＩは、互いに異なる中心波長λを有するＮ個のバンド（観測波長帯）Ｂ_１，Ｂ_２，…，Ｂ_Ｎにそれぞれ対応するＮ個のバンド画像ＩＢ_１，ＩＢ_２，…，ＩＢ_Ｎで構成されている。本実施の形態のバンド数Ｎは３以上の整数であるが、これに限定されずにバンド数Ｎが２であってもよい。たとえば、マルチスペクトルセンサの場合、バンド数Ｎは４程度であり、ハイパースペクトルセンサの場合、バンド数Ｎは数十〜数百程度である。これらバンド画像ＩＢ_１〜ＩＢ_Ｎは、すべて同一の被撮像領域を表している。また、バンド画像ＩＢ_１〜ＩＢ_Ｎの各画素は、２次元座標値ｘ，ｙの組で特定することができる。本実施の形態では、バンド画像ＩＢ_１〜ＩＢ_Ｎの各画素値は、放射輝度の観測値（単位：Ｗ／ｓｒ／ｍ^２）を示すものとする。

図３は、光学センサを搭載するプラットフォームが観測衛星３０である場合に、この光学センサによって惑星の地表が観測される様子を模式的に示す図である。図３に示されるように、観測衛星３０の光学センサは、地表の被撮像領域の方向から入射した電磁波を観測することにより、バンドＢ_ｉ（番号ｉは１〜Ｎのいずれか）の放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ，ｙ）（単位：Ｗ／ｓｒ／ｍ^２）を検出することができる。本実施の形態では、放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ，ｙ）は、近似的に次式（１）で表現されるものとする。

Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ，ｙ）
＝τ_ｉ×Ｌ_ｉ（ｘ，ｙ）＋α_ｉ（ｘ，ｙ）×Ｌ_{ｓｃａｔｔ０，ｉ} （１）

ここで、τ_ｉ×Ｌ_ｉ（ｘ，ｙ）は、地表から放射された後に大気を透過して光学センサに到達した電磁波の観測成分である。以下、この観測成分を「大気透過量」と呼ぶ。Ｌ_ｉ（ｘ，ｙ）は、地表の放射輝度であり、τ_ｉは、バンドＢ_ｉに依存する大気透過率である。

また、上式（１）におけるα_ｉ（ｘ，ｙ）×Ｌ_{ｓｃａｔｔ０，ｉ}は、大気で散乱された後に光学センサに入射した太陽光の観測成分である。以下、この観測成分を「大気散乱量」と呼ぶ。また、Ｌ_{ｓｃａｔｔ０，ｉ}は、バンドＢ_ｉと大気条件とに依存する成分である。以下、この成分を「散乱特性値」と呼ぶ。更に、α_ｉ（ｘ，ｙ）は、バンドＢ_ｉと被撮像領域の座標値ｘ，ｙとに依存する成分である。以下、この成分を「散乱度」と呼ぶ。

光線の大気散乱の原因としては、レイリー散乱及びミー散乱が支配的である。大気を伝搬する光線の波長と大気中の分子及び浮遊微粒子の構成とによって、大気散乱の大きさが決定される。このような大気散乱が生ずる実環境をシミュレーション計算用にモデル化することができる。本実施の形態では、観測対象となる惑星の大気散乱モデルに基づくシミュレーション計算によって、散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}が得られる。たとえば、公知のＭＯＤＴＲＡＮ（ＭＯＤｅｒａｔｅ ｒｅｓｏｌｕｔｉｏｎ ａｔｍｏｓｐｈｅｒｉｃ ＴＲＡＮｓｍｉｓｓｉｏｎ）と呼ばれる光波大気伝搬計算シミュレータを用いて散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}を計算することが可能である。一方、大気透過率τ_ｉとしては、大気透過モデルに基づくシミュレーション計算により得た値、もしくは、様々な環境下で測定された値を使用することができる。

図１を参照すると、画像処理装置１は、散乱特性データベース２５Ａ及び透過特性データベース２５Ｂが格納されたデータ記憶部２５を備えている。散乱特性データベース２５Ａには、各バンドＢ_ｉに対応する散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}が記憶されており、透過特性データベース２５Ｂには、各バンドＢ_ｉに対応する大気透過率τ_ｉが記憶されている。データ記憶部２５は、不揮発性メモリを用いて構成すればよい。

画像処理装置１は、更に、バンド画像ＩＢ_１〜ＩＢ_Ｎのうちの少なくとも１つのバンド画像の画素値のヒストグラムを算出するヒストグラム算出部１１と、当該算出されたヒストグラムに基づいて地表の被撮像領域における暗部領域を検出する暗部検出部１２と、散乱特性データベース２５Ａを用いて当該暗部領域に対応する大気の散乱度を推定する散乱度推定部１３と、当該推定された散乱度が予め定められた条件に照らして適正であるか否かを判定する散乱度判定部１４と、適正である判定された散乱度を用いて、被撮像領域における当該暗部領域以外の他の領域に対応する散乱度を補間する散乱度補間部１５と、この散乱度補間部１５で生成された散乱度分布及び透過特性データベース２５Ｂを用いてバンド画像ＩＢ_１〜ＩＢ_Ｎの画素値を補正する補正画像算出部１６と、補正されたバンド画像を出力するデータ出力部２１とを備えている。

以下、本実施の形態の画像処理装置１の構成例についてより詳細に説明する。図１に示されるヒストグラム算出部１１は、バンド画像ＩＢ_１〜ＩＢ_Ｎの中から暗部領域検出用のバンド画像を選択し、当該選択されたバンド画像の画素値のヒストグラムを算出する。暗部検出部１２は、このヒストグラムの頻度分布を正規分布などの確率分布で近似し、この確率分布を用いて暗部領域を検出（具体的には、暗部領域の座標値ｘ，ｙを検出）することができる。

霧や靄などの気象現象が発生していない領域のバンド画像からは、たとえば、図４のグラフに示すようなヒストグラムＨ_Ｒ，Ｈ_Ｇ，Ｈ_Ｂを得ることができる。このグラフでは、横軸が画素値に対応し、縦軸が各画素値の頻度に対応する。また、ヒストグラムＨ_Ｒは、赤色波長域のバンド画像の頻度分布を、ヒストグラムＨ_Ｇは、緑色波長域のバンド画像の頻度分布を、ヒストグラムＨ_Ｂは、青色波長域のバンド画像の頻度分布をそれぞれ示している。海や影などの低反射な被写体のバンド画像は、図５に示すように正規分布に近い分布のヒストグラムＨｘを持つと仮定することができる。

暗部検出部１２は、図５のヒストグラムＨｘの極大値のうち画素値が最も小さい極大値を中心とする頻度分布に正規分布をフィッティングすることにより正規分布Ｎｄを定めるパラメータを推定し、当該頻度分布を正規分布Ｎｄで近似することができる。フィッティングの方法としては、ガウス・ニュートン（Ｇａｕｓｓ−Ｎｅｗｔｏｎ）法またはマルカール（Ｍａｒｑｕａｒｄｔ）法などの公知の非線形最小自乗法を使用することができる。ここで、正規分布Ｎｄを定めるパラメータは、平均画素値μ及び分散σ^２（σは標準偏差）である。そして、暗部検出部１２は、この正規分布Ｎｄの裾部分の特定レベルＤＮ０以下の画素値を有する画素を暗部領域の画素とすることができる。特定レベルＤＮ０は、たとえば、正規分布Ｎｄの上側パーセント点（たとえば、上側５％点）に設定されればよい。なお、正規分布以外の他の確率分布を使用してフィッティングを実行してもよい。

一方、光学センサが、霧や靄などの気象現象が発生している領域を観測する場合には、たとえば、図６に示されるような赤色波長域、緑色波長域及び青色波長域のヒストグラムＨ_Ｒ，Ｈ_Ｇ，Ｈ_Ｂを有するバンド画像が得られる。この場合、ヒストグラムＨ_Ｒ，Ｈ_Ｇ，Ｈ_Ｂから暗部領域検出用の正規分布を得ることが難しいことがある。その理由は、霧や靄などの発生領域から伝搬した散乱光の強度が、地表面から伝搬した反射光の強度と同程度またはこれ以上となり、当該散乱光の観測成分と当該反射光の観測成分とが画素値に混在するためである。そこで、このような場合は、最も散乱の影響を受けない長波長側の赤外線バンドのバンド画像を用いて暗部領域の検出を行うことが望ましい。図６には、ＮＩＲ（近赤外線）バンドのバンド画像のヒストグラムＨ_ＮＩＲの例が示されている。なお、ＮＩＲバンドに限らず、水蒸気（Ｈ_２Ｏ）や二酸化炭素（ＣＯ_２）による吸収が少なく、霧、靄及び煙などを透過する赤外線のバンド画像を使用して暗部領域を検出してもよい。

散乱度推定部１３は、暗部検出部１２で検出された暗部領域の画素値Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ_ｄ，ｙ_ｄ）に基づき、散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}を用いて当該暗部領域に対応する大気の散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）を推定することができる。ここで、（ｘ_ｄ，ｙ_ｄ）は、当該暗部領域を特定する座標値であり、添え字「ｄ」は、当該暗部領域を指す記号である。暗部領域は、霧や靄などの発生領域に相当し、大気透過率τ_ｉが著しく低下した領域であるとみなされる。よって、上式（１）において、暗部領域については、大気透過量τ_ｉ×Ｌ_ｉ（ｘ_ｄ，ｙ_ｄ）を無視することができる。すなわち、暗部領域については、バンドＢ_ｉごとに、次式（１Ｄ）が近似的に成立する。

Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ_ｄ，ｙ_ｄ）＝α_ｉ（ｘ_ｄ，ｙ_ｄ）×Ｌ_{ｓｃａｔｔ０，ｉ} （１Ｄ）

散乱度推定部１３は、暗部領域の座標値（ｘ_ｄ，ｙ_ｄ）の各々について、バンド（波長）を説明変数として暗部領域の画素値Ｌ_{ｓｅｎｓｏｒ，１}（ｘ_ｄ，ｙ_ｄ）〜Ｌ_{ｓｅｎｓｏｒ，Ｎ}（ｘ_ｄ，ｙ_ｄ）を回帰分析することができる。この回帰分析により、暗部領域の座標値（ｘ_ｄ，ｙ_ｄ）の各々について回帰曲線Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）を得ることができる。

次に、散乱度推定部１３は、散乱特性データベース２５Ａから、各バンドＢ_ｉに対応する散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}を取得する。そして、散乱度推定部１３は、次式（２）に示されるように、各バンドＢ_ｉごとに、当該回帰曲線の値Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）を散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}で除算することにより、当該暗部領域に対応する大気の散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）を算出する。

α_ｉ（ｘ_ｄ，ｙ_ｄ）＝Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）／Ｌ_{ｓｃａｔｔ０，ｉ} （２）

図７は、暗部領域の特定の座標値（ｘ_ｄ，ｙ_ｄ）についてのバンド（観測波長）に対する大気散乱放射輝度の観測値（すなわち画素値）Ｌ_{ｓｅｎｓｏｒ，ｉ}（ｘ_ｄ，ｙ_ｄ）の例を示すグラフである。図７には、暗部領域におけるＮＩＲバンド画像の観測値Ｌ_ＮＩＲ（＝Ｌ_{ｓｅｎｓｏｒ，４}（ｘ_ｄ，ｙ_ｄ））、赤色バンド画像の観測値Ｌ_Ｒ（＝Ｌ_{ｓｅｎｓｏｒ，３}（ｘ_ｄ，ｙ_ｄ））、緑色バンド画像の観測値Ｌ_Ｇ（＝Ｌ_{ｓｅｎｓｏｒ，２}（ｘ_ｄ，ｙ_ｄ））及び青色バンド画像の観測値Ｌ_Ｂ（＝Ｌ_{ｓｅｎｓｏｒ，１}（ｘ_ｄ，ｙ_ｄ））が示されている。散乱度推定部１３は、これら観測値Ｌ_ＮＩＲ，Ｌ_Ｒ，Ｌ_Ｇ，Ｌ_Ｂを回帰分析して回帰曲線Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）を算出することができる。たとえば、最小自乗法を用いて、多項式関数で表現される回帰曲線Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）を算出すればよい。このとき、散乱度推定部１３は、上式（２）に従い、暗部領域の座標値（ｘ_ｄ，ｙ_ｄ）に対応する大気の散乱度α_１（ｘ_ｄ，ｙ_ｄ），α_２（ｘ_ｄ，ｙ_ｄ），α_３（ｘ_ｄ，ｙ_ｄ），α_４（ｘ_ｄ，ｙ_ｄ）を算出することができる。このように回帰曲線Ｆ（Ｂ_ｉ；ｘ_ｄ，ｙ_ｄ）を用いることにより、暗部領域に対応する大気の散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）を高い精度で算出することができる。

上述したように、散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}は、大気散乱モデルに基づき、ＭＯＤＴＲＡＮなどの光波大気伝搬計算シミュレータを用いて算出された値である。このため、本実施の形態では、暗部領域については、大気散乱モデルに基づいた大気散乱量α_ｉ（ｘ_ｄ，ｙ_ｄ）×Ｌ_{ｓｃａｔｔ０，ｉ}を算出することができる。これにより、大気散乱量α_ｉ（ｘ_ｄ，ｙ_ｄ）×Ｌ_{ｓｃａｔｔ０，ｉ}の算出に当たり、観測対象物（被写体）の分光反射特性の影響や雑音の影響を減らすことができる。

なお、大気散乱モデルに基づく散乱特性値Ｌ_{ｓｃａｔｔ０，ｉ}に代えて、教師データから得られた値を使用してもよい。

次に、散乱度判定部１４は、散乱度推定部１３で推定された散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）が所定条件に照らして適正であるか否かを判定する機能を有する。言い換えれば、散乱度判定部１４は、異常値を示す散乱度を除外する機能を有している。たとえば、散乱度判定部１４は、散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）が負の値である場合には、当該散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）が適正ではないと判定することができる。また、散乱度判定部１４は、大気散乱量α_ｉ（ｘ_ｄ，ｙ_ｄ）×Ｌ_{ｓｃａｔｔ０，ｉ}と、実際に観測された画素値とを互いに比較し、両者の差異が大きい場合（たとえば、両者の差異が閾値を超える場合）には、当該散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）は、大気散乱モデルから乖離した値であり、適正ではないと判定してもよい。あるいは、散乱度判定部１４は、被撮像領域全体に対応する散乱度分布に対して当該散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）が著しく異なる値を示す場合（たとえば、散乱度の平均値μについてμ±３σ以上の値を示す場合）に、当該散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）は適正ではないと判定することもできる。ここで、σは、散乱度分布の標準偏差である。

次に、散乱度補間部１５は、散乱度判定部１４で適正であると判定された散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）のみに基づいて、バンドＢ_ｉごとに、被撮像領域における暗部領域以外の他の領域に対応する散乱度を補間することにより、被撮像領域全体に対応する大気の散乱度分布を取得する。このとき、最近接（Ｎｅａｒｅｓｔ Ｎｅｉｇｈｂｏｒ）法、バイリニア（Ｂｉｌｉｎｅａｒ）法もしくはバイキュービック（Ｂｉｃｕｂｉｃ）法、またはこれらの組み合わせといった公知の画素補間法を使用して散乱度α_ｉ（ｘ，ｙ）を補間することができる。ここで、最近接法は、着目画素に最も近い画素位置の画素値をそのまま当該着目画素の画素値として使用する補間法である。また、バイリニア法は、着目画素の周辺４画素の画素値に基づき、１次多項式を用いて当該着目画素の画素値を算出する補間法である。そして、バイキュービック法は、着目画素の周辺１６画素の画素値に基づき、３次多項式を用いて当該着目画素の画素値を算出する補間法である。バイリニア法やバイキュービック法（特に、バイキュービック法）を使用すれば、連続的且つ滑らかな散乱度分布を推定することができる。

補正画像算出部１６は、散乱度補間部１５で得られた散乱度分布及び透過特性データベース２５Ｂを用いてバンド画像ＩＢ_１〜ＩＢ_Ｎの画素値を補正することにより、バンドＢ_１〜Ｂ_Ｎの補正画像を算出する。たとえば、次式（３）を用いて、バンドＢ_ｉの補正画像の画素値Ｄ_ｉ（ｘ，ｙ）を算出することができる。

Ｄ_ｉ（ｘ，ｙ）＝(Ｌ_ｉ（ｘ，ｙ）−Ｓ_ｉ（ｘ，ｙ）)／(τ_ｉ−Ｋ_ｉ（ｘ，ｙ）) （３）

ここで、Ｓ_ｉ（ｘ，ｙ）は、次式で定義される大気散乱量である。
Ｓ_ｉ（ｘ，ｙ）＝α_ｉ（ｘ，ｙ）×Ｌ_{ｓｃａｔｔ０，ｉ}

上式（３）の右辺の分母には、バンド画像のコントラストを補正するための調整項Ｋ_ｉ（ｘ，ｙ）が導入されている。この調整項Ｋ_ｉ（ｘ，ｙ）は、たとえば、次式（４）で近似的に表すことができる。

Ｋ_ｉ（ｘ，ｙ）＝Ｓ_ｉ（ｘ，ｙ）×β （４）

ここで、βは係数である。

大気透過率τ_ｉが増加すれば、その大気が光線を散乱させる量は低下する。逆に、霧や靄などの発生により大気透過率τ_ｉが低下すれば、その大気が光線を散乱させる量が増大して観測画像（バンド画像）のコントラストを低下させる状況が生ずると考えられる。そこで、上式（３）では、そのような状況を改善するために、大気散乱量Ｓ_ｉ（ｘ，ｙ）に比例する項Ｋ_ｉ（ｘ，ｙ）が導入されている。これにより、バンド画像ＩＢ_１〜ＩＢ_Ｎのコントラストを適正に補正することができる。

そして、データ出力部２１は、補正画像算出部１６で生成された補正画像を外部に出力する。このデータ出力部２１は、外部機器に補正画像のデータを転送する機能を有していてもよいし、あるいは、通信回線またはコンピュータネットワークを介して外部の通信機器（図示せず）へ補正画像のデータを送信する機能を有していてもよい。

次に、図８のフローチャートを参照しつつ、実施の形態１に係る画像処理の手順を説明する。

まず、画像入力部１０は、光学センサを搭載するプラットフォームと通信して、当該プラットフォームからマルチスペクトル画像ＭＳＩを取得する（ステップＳＴ１０）。ここで、画像入力部１０は、当該プラットフォームからマルチスペクトル画像ＭＳＩを取得する代わりに、マルチスペクトル画像ＭＳＩのデータが蓄積された地上のデータサーバにアクセスして、このデータサーバからマルチスペクトル画像ＭＳＩを取得してもよい。ヒストグラム算出部１１は、マルチスペクトル画像ＭＳＩを構成するバンド画像ＩＢ_１〜ＩＢ_Ｎのうちの少なくとも１つのバンド画像のヒストグラムを算出する（ステップＳＴ１１）。次いで、暗部検出部１２は、当該ヒストグラムに基づいて暗部領域を検出する（ステップＳＴ１２）。

次に、散乱度推定部１３は、散乱特性データベース２５Ａから散乱特性データを取得し（ステップＳＴ１３）、この散乱特性データを用いて当該暗部領域に対応する大気の散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）を推定する（ステップＳＴ１４）。散乱度判定部１４は、当該推定された散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）のうち適正な散乱度を抽出する（ステップＳＴ１５）。

次に、散乱度補間部１５は、当該適正な散乱度に基づいて、被撮像領域における当該暗部領域以外の他の領域に対応する散乱度を補間することにより、被撮像領域全体に対応する大気の散乱度分布を生成する（ステップＳＴ１６）。その後、補正画像算出部１６は、透過特性データベース２５Ｂから透過特性データを取得し（ステップＳＴ１７）、当該透過特性データ及び散乱度分布とを用いてバンド画像ＩＢ_１〜ＩＢ_Ｎの画素値を補正することにより、複数の補正画像を算出する（ステップＳＴ１８）。最終的に、データ出力部２１は、複数の補正画像を示す画像データを外部に出力する（ステップＳＴ１９）。

ところで、上記した画像処理装置１は、たとえば、ワークステーションまたはメインフレームなどの、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）内蔵のコンピュータで構成され得る。また、画像処理装置１の構成要素１０〜１６，２１の機能の全部または一部は、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）やＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）などの半導体集積回路で構成されてもよいし、あるいは、ＣＰＵを含むマイクロコンピュータの一種であるワンチップマイコンで構成されてもよい。

図９は、画像処理装置１のハードウェア構成例である情報処理装置１Ｓの概略構成を示す機能ブロック図である。図９に示されるように、この情報処理装置１Ｓは、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）４０Ｃを含むプロセッサ４０、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）４１、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）４２、大容量記録媒体４３及び入出力インタフェース４４を備えて構成されている。また、これらプロセッサ４０、ＲＡＭ４１、ＲＯＭ４２、大容量記録媒体４３及び入出力インタフェース４４は、バス４５を介して相互に接続されている。

ここで、入出力インタフェース４４は、図１に示した画像入力部１０及びデータ出力部２１に相当し、大容量記録媒体４３は、図１に示したデータ記憶部２５に相当する。入出力インタフェース４４は、外部機器（図示せず）と接続してデータの授受を行うためのデータ転送機能及び通信機能を有するものである。プロセッサ４０は、ＲＯＭ４２または大容量記録媒体４３からプログラムをロードし、当該プログラムに従って動作することによって画像処理装置１の機能を実現することができる。大容量記録媒体４３としては、たとえば、ハードディスク（磁気ディスク）、光ディスクまたはフラッシュメモリが挙げられる。

以上に説明したように実施の形態１の画像処理装置１では、散乱度推定部１３は、暗部領域に対応する大気の散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）を推定し、散乱度補間部１５は、当該散乱度α_ｉ（ｘ_ｄ，ｙ_ｄ）に基づいて被撮像領域全体に対応する散乱度分布を生成する。そして、補正画像算出部１６は、この散乱度分布を用いてバンド画像ＩＢ_１〜ＩＢ_Ｎを画素単位で補正することができる。したがって、霧や靄などによる太陽光の大気散乱の影響がバンド画像ＩＢ_１〜ＩＢ_Ｎ内に空間的に一様に現れない場合でも、各バンド画像を適正に補正することができる。これにより、より視認性の高いバンド画像を提供することが可能である。

実施の形態２．
次に、本発明に係る実施の形態２について説明する。図１０は、実施の形態２の画像処理装置２の概略構成を示す機能ブロック図である。図１０に示されるように、画像処理装置２は、上記実施の形態１の画像処理装置１と同様に、画像入力部１０、ヒストグラム算出部１１、暗部検出部１２、散乱度推定部１３、散乱度判定部１４、散乱度補間部１５及び補正画像算出部１６を備えている。これら構成要素１０〜１６の構成及び機能は、上述した通りである。

また、本実施の形態の画像処理装置２は、補正画像算出部１６で生成された補正画像を用いて、バンドＢ_ｉごとに被撮像領域の反射率マップρ_ｉ（ｘ，ｙ）を生成する反射率解析部１７と、この反射率マップρ_ｉ（ｘ，ｙ）を分析する分光分析部１８とを備えている。

更に、画像処理装置２は、データ記憶部２５に記憶された照度データベース２５Ｃを備える。照度データベース２５Ｃには、バンドＢ_ｉごとの太陽照度データが記憶されている。太陽照度データとしては、複数の撮像条件下で測定されたデータを使用することができる。あるいは、外部シミュレータを用いたシミュレーション計算により太陽照度データを算出してもよい。太陽照度は、撮像場所の大気粒子、エアロゾル及び太陽天頂角などの撮像条件によって変化する。そこで、たとえばＭＯＤＴＲＡＮを用いて複数の撮像条件下の太陽照度を計算し、これら太陽照度を示すデータを照度データベース２５Ｃに記録することができる。あるいは、外部シミュレータから照度データベース２５Ｃに太陽照度データが入力されてもよい。

反射率解析部１７は、次式（５）により、反射率マップρ_ｉ（ｘ，ｙ）を算出することができる。

ρ_ｉ（ｘ，ｙ）＝Ｄ_ｉ（ｘ，ｙ）／（Ｅ_ｉ／π） （５）

ここで、Ｅ_ｉは、バンドＢ_ｉに依存する太陽照度である。

分光分析部１８は、反射率解析部１０で算出された反射率マップρ_ｉ（ｘ，ｙ）をバンドＢ_ｉごとに分析する機能を持つ。分光分析部１８は、たとえば、正規化植生指標ＮＤＶＩ（Ｎｏｒｍａｌｉｚｅｄ Ｄｉｆｆｅｒｅｎｃｅ Ｖｅｇｅｔａｔｉｏｎ Ｉｎｄｅｘ）を用いて画像中の植生分布の把握や植生の分類を行ったり、正規化水指標（ＮＤＷＩ）を用いてバンド画像ＩＢ_１〜ＩＢ_Ｎ中の水域または非水域を検出したりすることができる。あるいは、分光分析部１８は、複数バンドＢ_１〜Ｂ_Ｎの反射率マップρ_１（ｘ，ｙ）〜ρ_Ｎ（ｘ，ｙ）を用いた主成分分析を行うことで、バンド画像ＩＢ_１〜ＩＢ_Ｎから水域、植生及び人工物などの画像データを抽出することもできる。

データ出力部２２は、分光分析部１８で生成された分析データを外部に出力する。このデータ出力部２２は、外部機器に分析データを転送する機能を有していてもよいし、あるいは、通信回線またはコンピュータネットワークを介して外部の通信機器（図示せず）へ分析データを送信する機能を有していてもよい。なお、データ出力部２２は、分析データに限らず、補正画像算出部１６で算出された補正画像群のデータを外部に出力するように構成されてもよい。

次に、図１１のフローチャートを参照しつつ、実施の形態２に係る画像処理の手順を説明する。図１１に示されるステップＳＴ１０〜ＳＴ１８の手順は、上記実施の形態１のステップＳＴ１０〜ＳＴ１８（図８）の手順と同じである。

ステップＳＴ１８の実行後、反射率解析部１７は、照度データベース２５Ｃから太陽照度データを取得し（ステップＳＴ２０）、この太陽照度データと補正画像の画素値とを用いて反射率マップρ_ｉ（ｘ，ｙ）を作成する（ステップＳＴ２１）。その後、分光分析部１８は、反射率マップρ_ｉ（ｘ，ｙ）を分析し、その分析結果を示す分析データをデータ出力部２２に与える（ステップＳＴ２２）。データ出力部２２は、分析データを外部に出力する（ステップＳＴ２３）。

なお、上記した画像処理装置２は、たとえば、ワークステーションまたはメインフレームなどの、ＣＰＵ内蔵のコンピュータで構成され得る。また、画像処理装置２の構成要素１０〜１８，２２の機能の全部または一部は、ＦＰＧＡやＡＳＩＣなどの半導体集積回路で構成されてもよいし、あるいは、ＣＰＵを含むマイクロコンピュータの一種であるワンチップマイコンで構成されてもよい。図９に示した情報処理装置１Ｓのハードウェア構成を用いて画像処理装置２を構成することもできる。

以上に説明したように実施の形態２の画像処理装置２は、補正画像算出部１６で得られた補正画像を用いて、地表における被撮像領域の反射率を解析することができるので、太陽の日照条件によらずに地表面の分光分析を行うことができる。

実施の形態３．
次に、本発明に係る実施の形態３について説明する。図１２は、実施の形態３の画像処理装置３の概略構成を示す機能ブロック図である。図１２に示されるように、画像処理装置３は、上記実施の形態１，２の画像処理装置１，２と同様に、画像入力部１０、ヒストグラム算出部１１、暗部検出部１２、散乱度推定部１３、散乱度判定部１４、散乱度補間部１５、補正画像算出部１６、反射率解析部１７及び分光分析部１８を備えている。これら構成要素１０〜１８の構成及び機能は、上述した通りである。

図１２に示されるように、画像処理装置３は、更に、分光分析部１８で生成された分析データに基づいて被撮像領域を分類する地表分類部１９と、その分類結果を示す分類データを外部に出力するデータ出力部２３とを備えている。地表分類部１９は、たとえば、ＮＤＶＩがある閾値以上の値を示す領域を植生領域に分類し、ＮＤＷＩが正の値を示す領域を水域に分類することができる。

データ出力部２３は、地表分類部１９で生成された分類データを外部に出力する。このデータ出力部２３は、外部機器に分類データを転送する機能を有していてもよいし、あるいは、通信回線またはコンピュータネットワークを介して外部の通信機器（図示せず）へ分類データを送信する機能を有していてもよい。なお、データ出力部２３は、分類データに限らず、補正画像算出部１６で算出された補正画像群のデータと、分光分析部１８で生成された分析データとを外部に出力するように構成されてもよい。

次に、図１３のフローチャートを参照しつつ、実施の形態３に係る画像処理の手順を説明する。図１３に示されるステップＳＴ１０〜ＳＴ１８の手順は、上記実施の形態１のステップＳＴ１０〜ＳＴ１８（図８）の手順と同じであり、図１３に示されるステップＳＴ２０〜ＳＴ２２の手順は、上記実施の形態２のステップＳＴ２０〜ＳＴ２２（図１１）の手順と同じである。ステップＳＴ２２の実行後、地表分類部１９は、分光分析部１８で生成された分析データに基づいて地表の被撮像領域を分類し、その分類結果を示す分類データをデータ出力部２３に与える（ステップＳＴ２４）。データ出力部２３は、分類データを外部に出力する（ステップＳＴ２５）。

なお、上記した画像処理装置３は、たとえば、ワークステーションまたはメインフレームなどの、ＣＰＵ内蔵のコンピュータで構成され得る。また、画像処理装置３の構成要素１０〜１９，２３の機能の全部または一部は、ＦＰＧＡやＡＳＩＣなどの半導体集積回路で構成されてもよいし、あるいは、ＣＰＵを含むマイクロコンピュータの一種であるワンチップマイコンで構成されてもよい。図９に示した情報処理装置１Ｓのハードウェア構成を用いて画像処理装置３を構成することもできる。

以上に説明したように実施の形態３の画像処理装置３は、分光分析部１８で生成された分析データに基づいて地表面の被覆分類を正確に行うことができる。

以上、図面を参照して本発明に係る種々の実施の形態１〜３について述べたが、これら実施の形態は本発明の例示であり、これら実施の形態以外の様々な形態を採用することもできる。本発明の範囲内において、上記実施の形態１〜３の自由な組み合わせ、各実施の形態の任意の構成要素の変形、または各実施の形態の任意の構成要素の省略が可能である。

ＩＢ_１〜ＩＢ_Ｎ バンド画像、１，２，３ 画像処理装置、１Ｓ 情報処理装置、１０ 画像入力部、１１ ヒストグラム算出部、１２ 暗部検出部、１３ 散乱度推定部、１４ 散乱度判定部、１５ 散乱度補間部、１６ 補正画像算出部、１７ 反射率解析部、１８ 分光分析部、１９ 地表分類部、２１〜２３ データ出力部、２５ データ記憶部、２５Ａ 散乱特性データベース、２５Ｂ 透過特性データベース、２５Ｃ 照度データベース、３０ 観測衛星、４０ プロセッサ、４０Ｃ ＣＰＵ、４１ ＲＡＭ、４２ ＲＯＭ、４３ 大容量記録媒体、４４ 入出力インタフェース、４５ バス。

Claims (12)

1. 大気を有する惑星の地表の被撮像領域を観測して複数の観測波長帯にそれぞれ対応する複数のバンド画像を生成する光学センサから、当該複数のバンド画像の供給を受ける画像入力部と、
   前記複数のバンド画像のうちの少なくとも１つのバンド画像の画素値に基づいて、前記被撮像領域における暗部領域を検出する暗部検出部と、
   前記観測波長帯に依存する散乱特性データが記憶されているデータ記憶部と、
   前記複数のバンド画像における当該暗部領域の画素値に基づき、前記散乱特性データを用いて当該暗部領域に対応する大気の散乱度を推定する散乱度推定部と、
   当該推定された散乱度に基づいて、前記被撮像領域における当該暗部領域以外の他の領域に対応する大気の散乱度を補間することにより、前記被撮像領域全体に対応する大気の散乱度分布を生成する散乱度補間部と、
   前記散乱度分布を用いて前記複数のバンド画像の画素値を補正することにより、前記複数の観測波長帯にそれぞれ対応する複数の補正画像を算出する補正画像算出部と
   を備えることを特徴とする画像処理装置。
2. 請求項１記載の画像処理装置であって、前記散乱度推定部は、前記複数のバンド画像における当該暗部領域の画素値を回帰分析して回帰曲線を算出し、当該回帰曲線の値を前記散乱特性データの値で除算することにより、当該暗部領域に対応する大気の散乱度を算出することを特徴とする画像処理装置。
3. 請求項１または請求項２記載の画像処理装置であって、前記散乱度推定部により推定された散乱度が適正であるか否かを判定する散乱度判定部を更に備え、
   前記散乱度補間部は、当該推定された散乱度のうち前記散乱度判定部により適正であると判定された散乱度に基づいて、前記他の領域に対応する大気の散乱度を補間することを特徴とする画像処理装置。
4. 請求項１から請求項３のうちのいずれか１項記載の画像処理装置であって、
   前記複数のバンド画像の各画素値は、前記大気で散乱された太陽光の観測成分である大気散乱量と、前記惑星の地表から放射された後に前記大気を透過した電磁波の観測成分である大気透過量とを含み、
   前記データ記憶部には、前記観測波長帯に依存する大気透過率を示す透過特性データが記憶されており、
   前記補正画像算出部は、前記散乱特性データの値に前記散乱度分布の散乱度を乗算することで大気散乱量を算出し、当該算出された大気散乱量及び前記透過特性データを用いて前記複数のバンド画像の画素値を補正する、
   ことを特徴とする画像処理装置。
5. 請求項１から請求項４のうちのいずれか１項記載の画像処理装置であって、前記暗部検出部は、前記複数のバンド画像のうち赤外線バンド画像の画素値に基づいて前記暗部領域を検出することを特徴とする画像処理装置。
6. 請求項１から請求項５のうちのいずれか１項記載の画像処理装置であって、当該少なくとも１つのバンド画像の画素値のヒストグラムを算出するヒストグラム算出部を更に備え、
   前記暗部検出部は、前記ヒストグラムに基づいて前記暗部領域を検出することを特徴とする画像処理装置。
7. 請求項６記載の画像処理装置であって、前記暗部検出部は、前記ヒストグラムの頻度分布を予め定められた確率分布で近似し、当該確率分布を用いて前記暗部領域を検出することを特徴とする画像処理装置。
8. 請求項１から請求項７のうちのいずれか１項記載の画像処理装置であって、
   前記複数の補正画像を用いて、前記観測波長帯ごとに前記被撮像領域の反射率マップを生成する反射率解析部と、
   前記反射率マップを分析する分光分析部と
   を更に備えることを特徴とする画像処理装置。
9. 請求項８記載の画像処理装置であって、前記分光分析部による分析結果に基づいて、前記被撮像領域を分類する地表分類部を更に備えることを特徴とする画像処理装置。
10. 大気を有する惑星の地表の被撮像領域を観測して複数の観測波長帯にそれぞれ対応する複数のバンド画像を生成する光学センサから当該複数のバンド画像の供給を受ける画像処理装置において実行される画像処理方法であって、
    前記複数のバンド画像のうちの少なくとも１つのバンド画像の画素値に基づいて前記被撮像領域における暗部領域を検出するステップと、
    前記観測波長帯に依存する散乱特性データが記憶されているデータ記憶部から、前記散乱特性データを取得するステップと、
    前記複数のバンド画像における当該暗部領域の画素値に基づき、前記散乱特性データを用いて当該暗部領域に対応する大気の散乱度を推定するステップと、
    当該推定された散乱度に基づいて、前記被撮像領域における当該暗部領域以外の他の領域に対応する大気の散乱度を補間することにより、前記被撮像領域全体に対応する大気の散乱度分布を生成するステップと、
    前記散乱度分布を用いて前記複数のバンド画像の画素値を補正することにより、前記複数の観測波長帯にそれぞれ対応する複数の補正画像を算出するステップと
    を備えることを特徴とする画像処理方法。
11. 請求項１０記載の画像処理方法であって、
    前記複数の補正画像を用いて、前記観測波長帯ごとに前記被撮像領域の反射率マップを生成するステップと、
    前記反射率マップを分析して分析データを出力するステップと
    を更に備えることを特徴とする画像処理方法。
12. 請求項１１記載の画像処理方法であって、前記分析データに基づいて、前記被撮像領域を分類するステップを更に備えることを特徴とする画像処理方法。

Images