JP6964834B2 - 画像処理装置および画像処理方法 - Google Patents

Description

本発明は、複数の観測波長帯の画像情報（以下、マルチスペクトル画像と記載する）を有した観測画像を処理する画像処理装置および画像処理方法に関する。

マルチスペクトルセンサまたはハイパースペクトルセンサと呼ばれる光学センサは、観測対象から放射された電磁波のうち、複数の観測波長帯の電磁波を検出することにより、観測波長帯ごとの画像情報を有した観測画像を生成することが可能である。例えば、この種の光学センサは、人工衛星または航空機に搭載されて、地表の観測に使用されている。

光学センサが高所からの地表の観測に使用される場合、光学センサによって検出される電磁波には、大気で散乱した太陽光の成分と、大気を透過してセンサに到達した電磁波の成分が含まれる。このため、光学センサによって生成された観測画像の分光特性は、地表で観測される分光特性とは異なるものとなる。従って、観測画像を用いて地表の分光特性を解析する場合、観測画像の成分から、大気の影響を受けた成分を除去する必要がある。

例えば、非特許文献１には、ＭＯＤＩＳ（Ｍｏｄｅｒａｔｅ Ｒｅｓｏｌｕｔｉｏｎ Ｉｍａｇｉｎｇ Ｓｐｅｃｔｒｏｒａｄｉｏｍｅｔｅｒ）と呼ばれる光学センサを用い、ＭＯＤＩＳによって生成された観測画像の成分から、大気の影響を受けた成分を除外する技術が記載されている。ＭＯＤＩＳとは、３６バンドの観測波長帯を有した光学センサであり、可視域から熱赤外域までの広い波長の範囲の電磁波を検出することが可能である。

非特許文献１に記載された従来の技術は、ＭＯＤＩＳによって観測可能な３６バンドの画像のうち、可視域から熱赤外域までの２０バンドの画像を用いてセンサの観測条件（エアロゾル、水蒸気量、オゾン量および雲領域など）を推定し、推定したデータを用いて大気の影響を低減させた可視域から赤外域までの７バンドの画像の分光反射率を算出する。大気の影響を低減させた分光反射率を用いることで、地表で観測される分光特性を持つ被写体（植生または土壌、人工物、水域など）の識別および分類が容易となり、あるいは、機械学習などを行うことで、観測条件に依存しない被写体の特徴量を把握することも可能である。

Ｅ． Ｆ． Ｖｅｒｍｏｔｅ ａｎｄ Ａ． Ｖｅｒｍｅｕｌｅｎ， " ＭＯＤＩＳ Ａｌｇｏｒｉｔｈｍ Ｔｅｃｈｎｉｃａｌ Ｂａｃｋｇｒｏｕｎｄ Ｄｏｃｕｍｅｎｔ ＡＴＭＯＳＰＨＥＲＩＣ ＣＯＲＲＥＣＴＩＯＮ ＡＬＧＯＲＩＴＨＭ：ＳＰＥＣＴＲＡＬ ＲＥＦＬＥＣＴＡＮＣＥＳ（ＭＯＤ０９） Ｖｅｒｓｉｏｎ ４．０ " ，Ａｐｒｉｌ １９９９， ＮＡＳＡ ｃｏｎｔａｃｔ ＮＡＳ５−９６０６２．

非特許文献１に記載された従来の技術は、地表面の分光反射率を推定するために用いる７バンドの画像に加えて、２０バンドの画像を大気の影響を推定するために用いるので、少なくとも２７バンドの画像を観測可能な光学センサが必要であった。ＭＯＤＩＳのように広い波長の範囲を多くのバンドで観測する光学センサは、観測条件ごとに適正なバンドを割り当てることはできる。しかしながら、観測可能なバンド数が多い分、光学センサを含む装置の規模が増大し、さらに、画像処理に使用されるバンド数が多いと、画像処理の負荷も高くなるという課題があった。

本発明は上記課題を解決するものであり、画像処理に使用される画像数および処理負荷を低減することができる画像処理装置および画像処理方法を得ることを目的とする。

本発明に係る画像処理装置は、地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数のマルチスペクトル画像を入力し、入力した複数のマルチスペクトル画像の各画素の画素値を放射輝度に変換する放射輝度変換部と、複数のマルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度をマルチスペクトル画像から検出する散乱特性検出部と、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部と、エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率のデータを抽出する大気特性抽出部と、大気特性抽出部によって抽出されたデータを用いて、被撮影領域の入射照度を算出する入射照度算出部と、マルチスペクトル画像の放射輝度値、大気特性抽出部によって抽出されたデータ、および入射照度算出部によって算出された入射照度を用いて、被撮影領域における地表面での分光反射率を算出する地表反射率算出部とを備える。

本発明によれば、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定し、推定したエアロゾルモデルから抽出された視程に対応する大気散乱放射輝度および大気透過率のデータと、マルチスペクトル画像の放射輝度値および入射照度を用いて、地表面での分光反射率を算出する。エアロゾルモデルを用いて大気散乱放射輝度および大気透過率が推定されるので、地表面の分光反射率の算出に使用するマルチスペクトル画像とは別に、大気の影響を推定するためのマルチスペクトル画像とその画像処理が不要である。このため、画像処理に使用される画像数および処理負荷を低減することができる。

実施の形態１に係る画像処理装置の構成を示す機能ブロック図である。 観測画像の構成を示す図である。 人工衛星に搭載された光学センサによって地表が観測される様子を模式的に示す図である。 実施の形態１に係る画像処理方法を示すフローチャートである。 観測画像から大気散乱放射輝度を検出する処理を示す概要図である。 図６Ａは、エアロゾルに対応する視程と大気散乱放射輝度との関係を示すグラフであり、図６Ｂは、視程に対応した大気散乱放射輝度の分光特性データを示すグラフであり、図６Ｃは、エアロゾルに対応した視程と大気透過率との関係を示すグラフである。 地表に入射される太陽光の成分を示す図である。 図８Ａは、実施の形態１に係る画像処理装置の機能を実現するハードウェア構成を示す機能ブロック図であり、図８Ｂは、実施の形態１に係る画像処理装置の機能を実現するソフトウェアを実行するハードウェア構成を示す機能ブロック図である。 実施の形態２に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態２に係る画像処理方法を示すフローチャートである。 実施の形態３に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態３に係る画像処理方法を示すフローチャートである。 地表と大気との間での太陽光の相互作用を示す図である。 実施の形態４に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態４に係る画像処理方法を示すフローチャートである。

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

図２は、観測画像ＭＳＩの構成を示す図である。観測画像ＭＳＩは、互いに異なる中心波長λを有するＮ個のバンド（観測波長帯）Ｂ_１，Ｂ_２，・・・，Ｂ_Ｎにそれぞれ対応するＮ枚のマルチスペクトル画像ＩＢ_１，ＩＢ_２，・・・，ＩＢ_Ｎから構成されている。なお、図２には、バンド数Ｎが３以上である場合を示したが、画像処理装置１は、これに限定されるものではなく、バンド数Ｎは２であってもよい。マルチスペクトル画像ＩＢ_１〜ＩＢ_Ｎは、全て同じ被撮影領域を表している。また、マルチスペクトル画像ＩＢ_１〜ＩＢ_Ｎの各画素は、２次元座標値ｘ，ｙの組みで特定することができる。

図３は、人工衛星１０に搭載された光学センサによって地表が観測される様子を模式的に示す図である。人工衛星１０に搭載された光学センサは、地表の被撮影領域の方向から入射した電磁波を観測することで、バンドＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}（単位；Ｗ／ｓｒ／ｍ^２／μｍ）を検出する。番号ｉは、１〜Ｎのいずれかである。また、放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}は、下記式（１）によって表現される。
Ｌ_{ｓｅｎｓｏｒ，ｉ}＝τ_ｉ×Ｌ_ｉ＋Ｌ_{ｓｃａｔｔ，ｉ} ・・・（１）

上記式（１）において、τ_ｉ×Ｌ_ｉは、地表から放射された後に大気を透過して光学センサに到達した電磁波の観測成分である。τ_ｉは、バンドＢ_ｉに依存する大気透過率であり、Ｌ_ｉは、地表の放射輝度である。上記式（１）におけるＬ_{ｓｃａｔｔ，ｉ}は、大気で散乱された光が光学センサに入射した成分である。

図１に示すように、画像処理装置１は、エアロゾルモデルデータベース１７Ａおよびスカイライトデータベース１７Ｂが格納されたデータ記憶部１７を備えている。エアロゾルモデルデータベース１７Ａには、大気中のエアロゾルに対応する視程と、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}および大気透過率τ_ｉの分光特性データとが対応付けて設定された複数のエアロゾルモデルが記憶されている。なお、エアロゾルは、大気中の浮遊微粒子である。視程は、大気中のエアロゾルによって、大気を伝搬する光が散乱または吸収されることで、その強度が一定割合（例えば、２％）の値に減少するまでの距離である。エアロゾルに対応する視程は、大気の特性を表す物理量である。

エアロゾルモデルにおいて、バンドＢ_ｉに依存した大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}（分光特性データ）およびバンドＢ_ｉに依存した大気透過率τ_ｉ（分光特性データ）のそれぞれの値が視程に対応している。エアロゾルモデルに設定される各データは、大気の影響を考慮した光波の伝搬環境をモデル化した大気散乱モデルおよび大気透過モデルに基づくシミュレーション計算によって得られる。例えば、ＭＯＤＴＲＡＮ（ＭＯＤｅｒａｔｅ ｒｅｓｏｌｕｔｉｏｎ ａｔｍｏｓｐｈｅｒｉｃ ＴＲＡＮｓｍｉｓｓｉｏｎ）と呼ばれる光波大気伝搬計算シミュレータを用いて算出することが可能である。なお、エアロゾルモデルには、様々な環境下で実測されたデータを設定することができる。

スカイライトデータベース１７Ｂは、エアロゾルモデルおよび視程とスカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}の関係が設定されたデータベースである。スカイライト（天空光）は、太陽からの直達光とは別に、大気中で散乱された後に地表に到達した太陽光、または、雲で反射された後に地表に到達した太陽光などである。スカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}は、例えば、ＭＯＤＴＲＡＮを用いたシミュレーション計算で得られた値、もしくは様々な環境下で実測された値を使用することができる。

さらに、画像処理装置１は、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６を備えている。放射輝度変換部１１は、光学センサからマルチスペクトル画像ＩＢ_ｉ（ｉ＝１〜Ｎ）を入力し、マルチスペクトル画像ＩＢ_ｉの各画素の画素値を放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}に変換する。散乱特性検出部１２は、マルチスペクトル画像ＩＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}の値に基づいて、マルチスペクトル画像ＩＢ_ｉから大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}を検出する。

エアロゾル・視程推定部１３は、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルの中から、散乱特性検出部１２によってマルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}に対応するエアロゾルモデルおよび視程を推定する。大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}および大気透過率τ_ｉの分光特性データを抽出する。

入射照度算出部１５は、被撮影領域の入射照度Ｅ_ｉを算出する。入射照度Ｅ_ｉは、太陽光の直達光成分の照度Ｅ_{ｄｉｒｅｃｔ，ｉ}とスカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}との和であり、バンドＢ_ｉに依存する。入射照度算出部１５は、スカイライトデータベース１７Ｂから、大気特性抽出部１４が大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}および大気透過率τ_ｉを抽出したときに用いたエアロゾルモデルおよび視程に対応するスカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}を取得して、入射照度Ｅ_ｉを算出する。地表反射率算出部１６は、放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}、大気透過率τ_ｉ、および入射照度Ｅ_ｉを用いて、被撮影領域における地表面での分光反射率を算出する。以下の説明では、地表面での分光反射率を、地表反射率ρ_ｉと記載する。

図４は、実施の形態１に係る画像処理方法を示すフローチャートであり、画像処理装置１による画像処理の手順を示している。まず、放射輝度変換部１１は、光学センサを搭載するプラットフォームと通信することで、プラットフォームに搭載された光学センサから観測画像ＭＳＩを受信する。また、放射輝度変換部１１は、プラットフォームに搭載された光学センサから観測画像ＭＳＩを入力する代わりに、観測画像ＭＳＩのデータが蓄積された地上のサーバにアクセスして、サーバから観測画像ＭＳＩを入力することもできる。

放射輝度変換部１１は、観測画像ＭＳＩを構成するマルチスペクトル画像ＩＢ_１〜ＩＢ_Ｎの各画素の画素値を、放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}に変換する（ステップＳＴ１）。例えば、マルチスペクトル画像ＩＢ_１〜ＩＢ_Ｎを構成している各画素の画素値は、ＤＮ（デジタルナンバー）と呼ばれる数値である。バンドＢ_ｉが可視域である場合に、ＤＮは、可視域を示すデジタル値ＤＮ（ｂａｎｄ）である。

マルチスペクトル画像ＩＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}は、下記式（２）を用いて算出することができる。α（ｂａｎｄ）は、光学センサにおける輝度値の感度に応じた放射輝度変換係数であり、β（ｂａｎｄ）は、光学センサにおける輝度値のオフセットに応じた放射輝度変換係数である。放射輝度変換部１１は、観測画像ＭＳＩに付随されたα（ｂａｎｄ）とβ（ｂａｎｄ）を取得し、取得したα（ｂａｎｄ）およびβ（ｂａｎｄ）を用いて、下記式（２）に従い、マルチスペクトル画像ＩＢ_ｉの各画素の画素値を放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}に変換する。
Ｌ_{ｓｅｎｓｏｒ，ｉ}＝α（ｂａｎｄ）×ＤＮ（ｂａｎｄ）＋β（ｂａｎｄ）・・・（２）

散乱特性検出部１２は、マルチスペクトル画像ＩＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}に基づいて、マルチスペクトル画像ＩＢ_ｉから、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}を検出する（ステップＳＴ２）。観察対象の地表に存在している被写体のうち、最も反射率が低い被写体を“低反射体”と呼ぶ。低反射体としては、例えば、海面などの水域または物体の陰影領域が挙げられる。散乱特性検出部１２は、マルチスペクトル画像ＩＢ_ｉの各画素の放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}から、低反射体に対応する領域の放射輝度を検出する。

図５は、観測画像ＭＳＩから、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}を検出する処理の概要を示す概要図である。散乱特性検出部１２は、マルチスペクトル画像ＩＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}をヒストグラム解析することで、図５中の右側のグラフを生成する。そして、散乱特性検出部１２は、マルチスペクトル画像ＩＢ_ｉを構成する画素のうち、低反射体に対応する領域の放射輝度を検出する。

低反射体が広範囲に一様かつ平坦な被写体（例えば、海面）である場合、地表面の分光反射率がほぼ０である。この場合、上記式（１）における大気透過量τ_ｉ×Ｌ_ｉが０になるので、マルチスペクトル画像ＩＢ_ｉから検出された放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}の値は、ほぼ大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}に一致する。この値は、マルチスペクトル画像から検出された推定値であるため、Ｌ’_{ｓｃａｔｔ，ｉ}と記載する。

次に、エアロゾル・視程推定部１３は、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルの中から、マルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ３）。なお、エアロゾルモデルには、様々な環境下における、エアロゾルに対応した視程と、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}および大気透過率τ_ｉの分光特性データとが対応付けて設定されている。

図６Ａは、エアロゾルに対応した視程と大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}との関係を示すグラフである。エアロゾルモデルには、例えば、図６Ａに示すデータがバンドＢ_ｉごとに設定されている。まず、エアロゾル・視程推定部１３は、散乱特性検出部１２によって検出された大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}を用いて、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルを参照する。そして、エアロゾル・視程推定部１３は、図６Ａに示す関係に基づいて、複数のエアロゾルモデルの中から、大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}に対応したエアロゾルモデルおよび大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}に対応した視程ｖ１を推定する。

次に、大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルから、推定された視程ｖ１に対応する大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}および大気透過率τ_ｉの分光特性データを抽出する（ステップＳＴ４）。図６Ｂは、視程に対応した大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}の分光特性データを示すグラフであり、バンドＢ１〜Ｂ_８のそれぞれにおける大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}を示している。図６Ｃは、エアロゾルに対応した視程ｖ１と大気透過率τ_ｉとの関係を示すグラフである。

大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルにおける図６Ｂに示す関係を参照することにより、視程ｖ１に対応した大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}の分光特性データを抽出する。ここで、データは離散的な視程の値に対して入力されるため、画像から推定した大気散乱放射輝度Ｌ’_{ｓｃａｔｔ，ｉ}とエアロゾルモデルデータベース１７Ａから抽出された大気散乱放射輝度（モデル値）Ｌ_{ｓｃａｔｔ，ｉ}を回帰分析し、エアロゾルモデルデータベース１７Ａから抽出された大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}に基づいた画像の大気散乱放射輝度値を算出する。例えば、大気散乱の大きさのみ変化するものとして、全バンドに固有の係数αを考慮したα×Ｌ_{ｓｃａｔｔ，ｉ}を画像の大気散乱放射輝度値とする。さらに、大気特性抽出部１４は、推定されたエアロゾルモデルにおける図６Ｃに示す関係を参照することにより、視程ｖ１に対応した大気透過率τ_ｉを抽出する。バンドＢ_ｉごとに推定された視程を用いて大気透過率τ_ｉの分光特性データが抽出される。

このように、画像処理装置１は、大気の影響の推定に使用されるデータ（大気散乱放射輝度および大気透過率）をエアロゾルモデルから推定するので、大気の影響を推定するためのマルチスペクトル画像が不要である。このため、観測画像ＭＳＩから地表反射率ρ_ｉを算出する画像処理に使用される画像数を低減することができ、これに伴い画像処理の回数が減るので、処理負荷も低減することができる。

次に、入射照度算出部１５は、被撮影領域の入射照度Ｅ_ｉを算出する（ステップＳＴ５）。図７は、地表に入射される太陽光の成分を示す図である。図７に示すように、地表には、太陽からの直達光に加え、スカイライトが到達する。入射照度算出部１５は、下記式（３）に従い、入射照度Ｅ_ｉを算出する。下記式（３）において、Ｅ_{ｄｉｒｅｃｔ，ｉ}は、バンドＢ_ｉに依存する直達光成分の照度であり、Ｅ_{ｓｋｙ，ｉ}は、バンドＢ_ｉに依存するスカイライト成分の照度である。なお、照度Ｅ_{ｄｉｒｅｃｔ，ｉ}は、大気上端（Ｔｏｐ Ｏｆ Ａｔｍｏｓｐｈｅｒｅ）の太陽照度Ｅ（ＴＯＡ）および太陽天頂角θから、下記式（４）を用いて算出することが可能である。
Ｅ_ｉ＝Ｅ_{ｄｉｒｅｃｔ，ｉ}＋Ｅ_{ｓｋｙ，ｉ} ・・・（３）
Ｅ_{ｄｉｒｅｃｔ，ｉ}＝τ_ｉ×Ｅ（ＴＯＡ）ｃｏｓθ ・・・（４）

入射照度算出部１５は、スカイライトデータベース１７Ｂを参照して、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルおよび視程に対応したスカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}を取得する。入射照度算出部１５は、観測画像ＭＳＩに付随した既定値である太陽照度Ｅ（ＴＯＡ）を取得し、太陽照度Ｅ（ＴＯＡ）、スカイライト成分の照度Ｅ_{ｓｋｙ，ｉ}、および大気特性抽出部１４によって抽出された大気透過率τ_ｉを用いて、上記式（３）および（４）に従い、入射照度Ｅ_ｉを算出する。

地表反射率算出部１６は、放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}、大気透過率τ_ｉおよび入射照度Ｅ_ｉを用いて、下記式（５）に従い、地表反射率ρ_ｉを算出する（ステップＳＴ６）。地表反射率ρ_ｉは、地表で観測される本来の分光反射率に近い値になるように大気の影響が低減されている。この地表反射率ρ_ｉを用いることで、観測画像ＭＳＩを、地表で観測される分光特性に合うように補正することが可能である。
ρ_ｉ＝（Ｌ_{ｓｅｎｓｏｒ，ｉ}−Ｌ_{ｓｃａｔｔ，ｉ}）／（τ_ｉ×Ｅ_ｉ／π） ・・・（５）

次に、画像処理装置１の機能を実現するハードウェア構成について説明する。
画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能は、処理回路によって実現される。すなわち、画像処理装置１は、図４のステップＳＴ１からステップＳＴ６までの処理を実行する処理回路を備える。処理回路は、専用のハードウェアであってもよいし、メモリに記憶されたプログラムを実行するＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）であってもよい。

図８Ａは、画像処理装置１の機能を実現するハードウェア構成を示す機能ブロック図であり、図８Ｂは、画像処理装置１の機能を実現するソフトウェアを実行するハードウェア構成を示す機能ブロック図である。図８Ａおよび図８Ｂにおいて、入出力インタフェース１００は、光学センサ（またはプラットフォーム）から画像処理装置１へ出力される観測画像ＭＳＩおよびこれに付随したデータを中継する。さらに、入出力インタフェース１００は、画像処理装置１から後段の装置（不図示）へ出力される地表反射率のデータを中継する。記憶装置１０１は、図１のデータ記憶部１７として機能する。

処理回路が、図８Ａに示す専用のハードウェアの処理回路１０２である場合、処理回路１０２は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能を別々の処理回路で実現してもよいし、これらの機能をまとめて１つの処理回路で実現してもよい。

処理回路が、図８Ｂに示すプロセッサ１０３である場合、画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能は、ソフトウェア、ファームウェアまたはソフトウェアとファームウェアの組み合わせによって実現される。ソフトウェアまたはファームウェアは、プログラムとして記述されてメモリ１０４に記憶される。

プロセッサ１０３は、メモリ１０４に記憶されたプログラムを読み出して実行することにより、画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能を実現する。例えば、画像処理装置１は、プロセッサ１０３によって実行されるときに、図４に示すフローチャートにおけるステップＳＴ１からステップＳＴ６の処理が結果的に実行されるプログラムを記憶するためのメモリ１０４を備える。これらのプログラムは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の手順または方法を、コンピュータに実行させる。メモリ１０４は、コンピュータを、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６として機能させるためのプログラムが記憶されたコンピュータ可読記憶媒体であってもよい。

メモリ１０４は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ−ＥＰＲＯＭ）などの不揮発性または揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤなどが該当する。

画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能の一部を専用ハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現してもよい。例えば、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３および大気特性抽出部１４は、専用のハードウェアである処理回路１０２によって機能を実現し、入射照度算出部１５および地表反射率算出部１６は、プロセッサ１０３がメモリ１０４に記憶されたプログラムを読み出し実行することで、機能を実現する。このように、処理回路は、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせにより、上記機能を実現することができる。

以上のように、実施の形態１に係る画像処理装置１は、観測画像ＭＳＩを構成するマルチスペクトル画像ＩＢ_ｉの画素値を放射輝度に変換する放射輝度変換部１１と、マルチスペクトル画像ＩＢ_ｉの放射輝度の値に基づいて、マルチスペクトル画像ＩＢ_ｉから大気散乱放射輝度を検出する散乱特性検出部１２と、複数のエアロゾルモデルの中から、マルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部１３と、推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率の分光特性データを抽出する大気特性抽出部１４と、大気透過率データを用いて、入射照度Ｅ_ｉを算出する入射照度算出部１５と、マルチスペクトル画像ＩＢ_ｉの放射輝度値、大気散乱放射輝度および大気透過率のデータおよび入射照度Ｅ_ｉを用いて、地表反射率ρ_ｉを算出する地表反射率算出部１６を備えている。エアロゾルモデルから大気散乱放射輝度および大気透過率が推定されるので、これらのデータを推定するためにマルチスペクトル画像および画像処理が不要である。これにより、画像処理に使用される画像数および処理負荷を低減することが可能である。

実施の形態２．
図９は、実施の形態２に係る画像処理装置１Ａの構成を示す機能ブロック図である。図９に示すように、画像処理装置１Ａは、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃおよびエアロゾルモデルデータベース１７Ｄが格納されたデータ記憶部１７を備えている。画像処理装置１Ａは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６および地表アルベド算出部１８を備えている。ここで、スカイライトデータベース１７Ｂ、放射輝度変換部１１、散乱特性検出部１２、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能および各構成は、実施の形態１で説明した通りである。

地表アルベドデータベース１７Ｃは、地表の位置情報とその位置における地表面のアルベドとが対応付けて設定されたデータベースである。地表面のアルベドは、太陽光が地球の地表面で反射される割合を示す指標値であり、反射率に相当する。以下、地表面のアルベドを“地表アルベド”と記載する。地表アルベドが高い場合、地表面で反射して大気中で散乱される太陽光の割合が増加するため、大気での散乱光が大きくなる。例えば、地表がコンクリートまたは砂漠領域であれば、地表アルベドが高く、これに伴って大気散乱光も大きくなる。一方、海面または陰影領域では、地表アルベドが低いため、大気散乱光も小さくなる。地表アルベドデータベース１７Ｃに設定される地表アルベドには、例えば、ＭＯＤＴＲＡＮを用いたシミュレーション計算で得られた値、もしくは、様々な環境下で実測された値を使用することができる。

エアロゾルモデルデータベース１７Ｄにおいて、複数のエアロゾルモデルのそれぞれが地表アルベドに対応付けられている。例えば、エアロゾルモデルには、エアロゾルおよび地表アルベドに対応する視程と、大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定されている。

地表アルベド算出部１８は、地表アルベドデータベース１７Ｃに格納されている地表アルベドを用いて、被撮影領域における地表アルベドを算出する。なお、観測画像ＭＳＩに撮影された地表で反射された光は、大気中で多重散乱されている。そこで、地表アルベド算出部１８は、観測画像ＭＳＩの被撮影領域に対応する全ての地表アルベドを地表アルベドデータベース１７Ｃから取得すると、取得した地表アルベドの平均値を算出する。地表アルベド算出部１８は、地表アルベドの平均値を、被撮影領域における地表アルベドとしてエアロゾル・視程推定部１３Ａに出力する。

エアロゾル・視程推定部１３Ａは、エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルのうち、地表アルベドおよび大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する。例えば、エアロゾル・視程推定部１３Ａは、地表アルベド算出部１８によって算出された地表アルベドを入力し、散乱特性検出部１２によってマルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度を入力して、入力した地表アルベドおよび大気散乱放射輝度に対応するエアロゾルモデルおよび視程を、エアロゾルモデルデータベース１７Ｄから推定する。

図１０は、実施の形態２に係る画像処理方法を示すフローチャートであって、画像処理装置１Ａによる画像処理の手順を示している。図１０のステップＳＴ１ａおよびステップＳＴ２ａは、図４のステップＳＴ１およびステップＳＴ２と同様の処理であり、ステップＳＴ５ａからステップＳＴ７ａまでの処理は、図４のステップＳＴ４からステップＳＴ６と同様の処理であるので、これらのステップの説明は省略する。

ステップＳＴ２ａの処理が終了した後、地表アルベド算出部１８は、被撮影領域の地表アルベドを算出する（ステップＳＴ３ａ）。地表アルベド算出部１８によって算出される地表アルベドは、例えば、観測画像ＭＳＩの被撮影領域の位置情報に基づいて、地表アルベドデータベース１７Ｃから取得された地表アルベドの平均値である。

エアロゾル・視程推定部１３Ａは、エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルのうち、地表アルベド算出部１８によって算出された地表アルベドおよび散乱特性検出部１２によって検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ４ａ）。エアロゾル・視程推定部１３Ａによって推定されたエアロゾルモデルおよび視程は、大気特性抽出部１４に出力される。

なお、画像処理装置１Ａにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６および地表アルベド算出部１８の各機能は、処理回路によって実現される。すなわち、画像処理装置１Ａは、図１０に示したステップＳＴ１ａからステップＳＴ７ａまでの処理を実行するための処理回路を備えている。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。

以上のように、実施の形態２に係る画像処理装置１Ａは、被撮影領域の地表アルベドを算出する地表アルベド算出部１８を備えている。エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルは、地表アルベドに対応付けられている。エアロゾル・視程推定部１３Ａは、地表アルベド算出部１８によって算出された地表アルベドおよびマルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する。被撮影領域の地表アルベドが考慮されたエアロゾルモデルおよび視程を推定することで、エアロゾルモデルおよび視程の推定精度が向上する。

実施の形態３．
図１１は、実施の形態３に係る画像処理装置１Ｂの構成を示す機能ブロック図である。図１１に示すように、画像処理装置１Ｂは、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃおよびエアロゾルモデルデータベース１７Ｄが格納されたデータ記憶部１７を備えている。さらに、画像処理装置１Ｂは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６Ａおよび地表アルベド算出部１８を備えている。また、地表反射率算出部１６Ａは、補正部１９を備えている。

なお、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃ、エアロゾルモデルデータベース１７Ｄ、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５および地表アルベド算出部１８の各機能および各構成は、実施の形態１および実施の形態２で説明した通りである。

地表反射率算出部１６Ａは、実施の形態１と同様の手順で、マルチスペクトル画像ＩＢ_ｉの放射輝度Ｌ_{ｓｅｎｓｏｒ，ｉ}、大気散乱放射輝度Ｌ_{ｓｃａｔｔ，ｉ}、大気透過率τ_ｉおよび入射照度Ｅ_ｉを用いて、被撮影領域の地表反射率ρ_ｉを算出する。補正部１９は、地表と大気との間での光の相互作用に応じた地表反射率ρ_ｉの変化分を補正する。

図１２は、実施の形態３に係る画像処理方法を示すフローチャートであって、画像処理装置１Ｂによる画像処理の手順を示している。図１２のステップＳＴ１ｂからステップＳＴ７ｂまでの処理は、図１０のステップＳＴ１ａからステップＳＴ７ａまでと同様の処理であるので、これらのステップの説明は省略する。

補正部１９は、地表と大気との間における太陽光の相互作用に応じた地表反射率ρ_ｉの変化分を補正する（ステップＳＴ８ｂ）。図１３は、地表と大気との間での太陽光の相互作用を示す図である。地表アルベドデータベース１７Ｃにおける地表アルベドが得られた環境と観測画像ＭＳＩが撮影された観測環境との差異に応じて大気散乱光の大きさは変化する。

例えば、地表反射率ρ_ｉの算出に使用された観測画像ＭＳＩの観測環境と、地表アルベドデータベース１７Ｃにおける地表アルベドの算出に使用された観測画像ＭＳＩの観測環境との間で、光学センサの視線方向または太陽高度に差異があると、大気散乱光の大きさまたは被写体のＢＲＤＦ（双方向散乱分布関数）が変化し得る。これにより、大気から地球への散乱光が変化し、地表を照らす照度が変化し得る。従って、入射照度Ｅ_ｉは、太陽光の直達光成分とスカイライト成分の組み合わせによってその分光特性が変化する。補正部１９は、地表反射率ρ_ｉを、地表へ入射されるスカイライト成分の変化の影響が低減された本来の値に近い地表反射率ρ_ｉ’に補正する。

補正部１９には、例えば、地表アルベドデータベース１７Ｃに格納された地表アルベドが閾値ρｔｈとして設定されている。補正部１９は、地表反射率ρ_ｉが閾値ρｔｈとの間で下記式（６）の関係にある場合、地表と大気との間で太陽光のスカイライト成分が相互作用して地表へ入射されるスカイライト成分が変化したと推定する。
ρ_ｉ≦ρｔｈ ・・・（６）

地表と大気との間での相互作用によるスカイライト成分の変化を、元のスカイライト成分からの割合として係数ｒ_ｉとした場合に、スカイライトの照度成分の変化分を含む入射照度Ｅ_ｉ’は、下記式（７）で表される。入射照度Ｅ_ｉ’を用いて算出される地表反射率ρ_ｉ’は、地表反射率ρ_ｉおよびその算出に使用された入射照度Ｅ_ｉを用いることで、下記式（８）に示す関係で表すことができる。補正部１９は、地表へ入射されるスカイライト成分が変化したと推定すると、下記式（７）に従って、地表反射率ρ_ｉにおけるスカイライト成分に応じた変化分を低減させた地表反射率ρ_ｉ’に補正する。これは、例えば、太陽高度が低く、建物などの影が多いような画像を想定している。逆にρ_ｉがρｔｈ以上である場合に、スカイライトの成分がより大きくなる（ｒ_ｉ＞１）ような補正を行ってもよい。これは、例えば、順光によって正反射成分が強く作用するような画像を想定している。
Ｅ_ｉ’＝Ｅ_{ｄｉｒｅｃｔ，ｉ}＋ｒ_ｉ×Ｅ_{ｓｋｙ，ｉ} ・・・（７）
ρ_ｉ’＝ρ_ｉ（Ｅ_ｉ／Ｅ_ｉ’） ・・・（８）

なお、画像処理装置１Ｂにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６Ａおよび地表アルベド算出部１８の各機能は、処理回路によって実現される。すなわち、画像処理装置１Ｂは、図１２に示したステップＳＴ１ｂからステップＳＴ８ｂまでの処理を実行するための処理回路を備える。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。

以上のように、実施の形態３に係る画像処理装置１Ｂは、地表と大気との間での光の相互作用に応じた地表反射率ρ_ｉの変化分を補正する補正部１９を備えている。地表と大気との間での相互作用によって地表へ入射されるスカイライト成分が変化しても、この影響を低減させた地表反射率ρ_ｉ’を得ることができる。

実施の形態４．
図１４は、実施の形態４に係る画像処理装置１Ｃの構成を示す機能ブロック図である。図１４に示すように、画像処理装置１Ｃは、スカイライトデータベース１７Ｂおよびエアロゾルモデルデータベース１７Ｅが格納されたデータ記憶部１７を備えている。さらに、画像処理装置１Ｃは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ｂ、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６を備えている。エアロゾル・視程推定部１３Ｂは、変換部１３Ｃを備える。なお、スカイライトデータベース１７Ｂ、放射輝度変換部１１、散乱特性検出部１２、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能および各構成は、実施の形態１で説明した通りである。

本発明の発明者による実験解析の結果、太陽天頂角θが小さい範囲では、太陽天頂角θが変化しても大気散乱放射輝度の絶対値の変化はあるものの、分光特性の変化は小さいことが判明した。すなわち、太陽天頂角θが小さい範囲では、太陽天頂角θに対する絶対値の補正係数を与えれば、エアロゾルモデルに設定されている大気散乱放射輝度を、太陽天頂角θが０度であるときの大気散乱放射輝度に置き換えても、置き換え前の大気散乱放射輝度との対応関係と大きく乖離しないことを意味する。なお、この特性は、光学センサを搭載する人工衛星の天頂角θ_ｓａｔが約４０度であるときには、太陽天頂角θが約５０度以下でこの関係を満たす。

実施の形態１で示したエアロゾルモデルデータベース１７Ａにおけるエアロゾルモデルには、大気散乱放射輝度および大気透過率の分光特性データが太陽天頂角θごとに設定されている。これに対して、エアロゾルモデルデータベース１７Ｅでは、太陽天頂角θが小さい範囲において、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に分光特性データを対応付ければよい。このため、エアロゾルモデルデータベース１７Ｅは、エアロゾルモデルデータベース１７Ａに比べて、太陽天頂角θが小さい範囲におけるデータの量が大幅に削減される。

エアロゾル・視程推定部１３Ｂには、例えば、太陽天頂角θに対する大気散乱放射輝度値の変化が、ルックアップテーブル（以下、ＬＵＴと記載する）として設定されている。太陽天頂角θが閾値以下である場合、変換部１３Ｃは、このＬＵＴを参照して、散乱特性検出部１２によって検出された大気散乱放射輝度を、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に変換する。エアロゾルモデル・視程推定部１３Ｂは、変換部１３Ｃによって変換された大気散乱放射輝度を、エアロゾルモデルおよび視程の推定に使用する。

図１５は、実施の形態４に係る画像処理方法を示すフローチャートであって、画像処理装置１Ｃによる画像処理の手順を示している。図１５のステップＳＴ１ｃ、ステップＳＴ２ｃおよびステップＳＴ６ｃからステップＳＴ８ｃの処理は、図４におけるステップＳＴ１、ステップＳＴ２およびステップＳＴ４からステップＳＴ６までと同様の処理であるので、これらのステップの説明は省略する。

ステップＳＴ２ｃの処理が終了すると、エアロゾル・視程推定部１３Ｂは、散乱特性検出部１２によって大気散乱放射輝度が検出された画像が観測されたときの太陽天頂角θが閾値以下であるか否かを判定する（ステップＳＴ３ｃ）。なお、太陽天頂角θは、観測画像ＭＳＩに付随している既知データである。また、閾値は、例えば６０度とする。

太陽天頂角θが閾値以下であった場合（ステップＳＴ３ｃ；ＹＥＳ）、変換部１３Ｃは、ＬＵＴを参照して、散乱特性検出部１２によって検出された大気散乱放射輝度を、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に変換する（ステップＳＴ４ｃ）。この後、エアロゾル・視程推定部１３Ｂは、エアロゾルモデルデータベース１７Ａから、変換部１３Ｃによって変換された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ５ｃ）。

一方、太陽天頂角θが閾値よりも大きかった場合（ステップＳＴ２ｃ；ＮＯ）、エアロゾル・視程推定部１３Ｂは、エアロゾルモデルデータベース１７Ａから、散乱特性検出部１２によってマルチスペクトル画像ＩＢ_ｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ５ｃ）。

なお、画像処理装置１Ｃにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ｂ、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能は、処理回路により実現される。すなわち、画像処理装置１Ｃは、図１５に示したステップＳＴ１ｃからステップＳＴ８ｃまでの処理を実行するための処理回路を備えている。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。

以上のように、実施の形態４に係る画像処理装置１Ｃにおいて、エアロゾル・視程推定部１３Ｂが、太陽天頂角が閾値以下である場合、マルチスペクトル画像から検出された大気散乱放射輝度が太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用する。エアロゾルモデルデータベース１７Ｅが、太陽天頂角θが小さい範囲において、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に分光特性データを対応付ければよい。このため、エアロゾルモデルデータベース１７Ｅにおける、太陽天頂角θが小さい範囲に対応するデータの量を大幅に削減することが可能である。

なお、本発明は上記実施の形態に限定されるものではなく、本発明の範囲内において、実施の形態のそれぞれの自由な組み合わせまたは実施の形態のそれぞれの任意の構成要素の変形もしくは実施の形態のそれぞれにおいて任意の構成要素の省略が可能である。

本発明に係る画像処理装置は、例えば、人工衛星または航空機に搭載された光学センサによって観測された観測画像の画像処理に利用可能である。

１，１Ａ，１Ｂ，１Ｃ 画像処理装置、１０ 人工衛星、１１ 放射輝度変換部、１２ 散乱特性検出部、１３，１３Ａ，１３Ｂ エアロゾル・視程推定部、１３Ｃ 変換部、１４ 大気特性抽出部、１５ 入射照度算出部、１６，１６Ａ 地表反射率算出部、１７ データ記憶部、１７Ａ，１７Ｄ，１７Ｅ エアロゾルモデルデータベース、１７Ｂ スカイライトデータベース、１７Ｃ 地表アルベドデータベース、１８ 地表アルベド算出部、１９ 補正部、１００ 入出力インタフェース、１０１ 記憶装置、１０２ 処理回路、１０３ プロセッサ、１０４ メモリ。

Claims (8)

1. 地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数の前記マルチスペクトル画像を入力し、入力した複数の前記マルチスペクトル画像の各画素の画素値を放射輝度に変換する放射輝度変換部と、
   複数の前記マルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度を前記マルチスペクトル画像から検出する散乱特性検出部と、
   エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部と、
   前記エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率のデータを抽出する大気特性抽出部と、
   前記大気特性抽出部によって抽出されたデータを用いて、前記被撮影領域の入射照度を算出する入射照度算出部と、
   前記マルチスペクトル画像の放射輝度値、前記大気特性抽出部によって抽出されたデータ、および前記入射照度算出部によって算出された入射照度を用いて、前記被撮影領域における地表面での分光反射率を算出する地表反射率算出部と、
   を備えたことを特徴とする画像処理装置。
2. 前記被撮影領域における地表面のアルベドを算出する地表アルベド算出部を備え、
   複数のエアロゾルモデルは、前記アルベドに対応付けられており、
   前記エアロゾル・視程推定部は、複数のエアロゾルモデルのうち、前記地表アルベド算出部によって算出された前記アルベドおよび前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定すること
   を特徴とする請求項１記載の画像処理装置。
3. 地表と大気との間での光の相互作用に応じた前記分光反射率の変化分を補正する補正部を備えたこと
   を特徴とする請求項２記載の画像処理装置。
4. 前記エアロゾル・視程推定部は、太陽天頂角が閾値以下である場合、前記マルチスペクトル画像から検出された大気散乱放射輝度が前記太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用すること
   を特徴とする請求項１から請求項３のいずれか１項記載の画像処理装置。
5. 放射輝度変換部が、地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数の前記マルチスペクトル画像を入力し、入力した複数の前記マルチスペクトル画像の各画素の画素値を放射輝度に変換するステップと、
   散乱特性検出部が、複数の前記マルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度を前記マルチスペクトル画像から検出するステップと、
   エアロゾル・視程推定部が、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するステップと、
   大気特性抽出部が、前記エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率の分光特性データを抽出するステップと、
   入射照度算出部が、前記大気特性抽出部によって抽出されたデータを用いて、前記被撮影領域の入射照度を算出するステップと、
   地表反射率算出部が、前記マルチスペクトル画像の放射輝度値、前記大気特性抽出部によって抽出されたデータ、および前記入射照度算出部によって算出された入射照度を用いて、前記被撮影領域における地表面での分光反射率を算出するステップと、
   を備えたことを特徴とする画像処理方法。
6. 地表アルベド算出部が、前記被撮影領域における地表面のアルベドを算出するステップを備え、
   複数のエアロゾルモデルは、前記アルベドに対応付けられており、
   前記エアロゾル・視程推定部は、複数のエアロゾルモデルのうち、前記地表アルベド算出部によって算出された前記アルベドおよび前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定すること
   を特徴とする請求項５記載の画像処理方法。
7. 補正部が、地表と大気との間での光の相互作用に応じた前記分光反射率の変化分を補正するステップを備えたこと
   を特徴とする請求項６記載の画像処理方法。
8. 前記エアロゾル・視程推定部が、太陽天頂角が閾値以下である場合、前記マルチスペクトル画像から検出された大気散乱放射輝度が前記太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用するステップを備えたこと
   を特徴とする請求項５から請求項７のいずれか１項記載の画像処理方法。

Images