JPWO2018193517A1 - 画像処理装置 - Google Patents

Abstract

散乱輝度算出部（３）により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部（７）を設け、補正部（１２）が、補正値算出部（７）により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する。

Description

この発明は、パンクロマチック画像における分光放射輝度を補正する画像処理装置に関するものである。

センサから取得した画像の一部分のコントラストが、大気中に含まれる靄又は霧などの影響を受けて低下していることがある。
コントラストの低下は、画像の視認性の低下を招くため、靄又は霧などの影響で低下しているコントラストを補正する画像処理装置が開発されている。

以下の特許文献１には、画像を構成する複数の画素のうち、任意の画素の信号レベルと、当該画素の周囲の領域に存在している画素の信号レベルとの差を計測し、その差が大きければ、コントラストが低下しているとみなして、当該画素の信号レベルを高める画像処理装置が開示されている。

特開２０１２−０２７５４７号公報

例えば、センサから取得した画像が海岸線を撮影した画像である場合、画像を構成する複数の画素のうち、任意の画素が海の領域に対応し、当該画素の周囲の領域に存在している画素が海岸の領域に対応していることがある。
このとき、海の反射率と海岸の反射率とは大きな差異があるため、大気中に含まれる靄又は霧などの影響を受けていない場合でも、任意の画素の信号レベルと、当該画素の周囲の領域に存在している画素の信号レベルとの差が大きくなることがある。
このような場合に、コントラストが低下しているとみなして、当該画素の信号レベルを高めてしまうと、却って、画像の視認性が低下する。
このため、センサから取得した画像が海岸線を撮影した画像である場合のほか、山間部と市街地の境界等のように、反射率が大きく異なる複数の領域を含んでいる画像である場合、靄又は霧などの影響を受けて低下しているコントラストを補正して、画像の視認性を高めることが困難であるという課題があった。

この発明は上記のような課題を解決するためになされたもので、センサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる画像処理装置を得ることを目的とする。

この発明に係る画像処理装置は、複数の波長帯の画像を含んでいるマルチスペクトル画像を複数の領域に分割する画像分割部と、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部により分割された各々の領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する散乱輝度算出部と、散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部とを設け、補正部が、補正値算出部により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するようにしたものである。

この発明によれば、散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部を設け、補正部が、補正値算出部により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、センサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる効果がある。

この発明の実施の形態１による画像処理装置を示す構成図である。 この発明の実施の形態１による画像処理装置を示すハードウェア構成図である。 画像処理装置がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。 画像処理装置がソフトウェア又はファームウェアなどで実現される場合の処理手順を示すフローチャートである。 画像分割部１により分割されたマルチスペクトル画像を示す説明図である。 画像分割部１により分割されたマルチスペクトル画像を示す説明図である。 暗部輝度算出部４により算出されたバンドｂにおける小領域Ｓ（ｉ，ｊ）のヒストグラムを示す説明図である。 散乱輝度算出処理部５によりモデル化された大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を示す説明図である。 バンドｂにおける小領域Ｓ（ｉ，ｊ）の大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）に対応する大気透過率τ（ｂ，ｉ，ｊ）を示す説明図である。 マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯を示す説明図である。 この発明の実施の形態２による画像処理装置を示す構成図である。 この発明の実施の形態２による画像処理装置を示すハードウェア構成図である。 反射輝度算出部４１によりモデル化された分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）を示す説明図である。 この発明の実施の形態３による画像処理装置を示す構成図である。 この発明の実施の形態３による画像処理装置を示すハードウェア構成図である。

以下、この発明をより詳細に説明するために、この発明を実施するための形態について、添付の図面に従って説明する。

実施の形態１．
図１は、この発明の実施の形態１による画像処理装置を示す構成図であり、図２は、この発明の実施の形態１による画像処理装置を示すハードウェア構成図である。
図１及び図２において、画像分割部１は、例えば図２に示す画像分割回路２１で実現される。
画像分割部１は、マルチスペクトルセンサから複数の波長帯の画像を含んでいるマルチスペクトル画像を取得し、マルチスペクトル画像を複数の小領域（領域）に分割する処理を実施する。
パラメータ記憶部２は、例えば図２に示すパラメータ記憶回路２２で実現される。
パラメータ記憶部２は、各々の波長帯の散乱特性を示す散乱特性データ、小領域の大気透過率を示す関数である散乱透過特性データ及びセンサの観測可能な波長帯を示す波長特性データを記憶している。

散乱輝度算出部３は、例えば図２に示す散乱輝度算出回路２３で実現され、暗部輝度算出部４及び散乱輝度算出処理部５を備えている。
散乱輝度算出部３は、パラメータ記憶部２により記憶されている散乱特性データを用いて、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する処理を実施する。

散乱輝度算出部３の暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域のヒストグラムをそれぞれ算出し、そのヒストグラムから、画像分割部１により分割された各々の小領域の中で相対的に分光放射輝度が小さい暗部領域の分光放射輝度をそれぞれ算出する処理を実施する。
散乱輝度算出部３の散乱輝度算出処理部５は、パラメータ記憶部２により記憶されている散乱特性データ及び暗部輝度算出部４により算出された暗部領域の分光放射輝度を用いて、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱放射輝度をそれぞれ算出する処理を実施する。

大気透過率算出部６は、例えば図２に示す透過率算出回路２４で実現される。
大気透過率算出部６は、パラメータ記憶部２により記憶されている散乱透過特性データを用いて、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱の大気透過率をそれぞれ算出する処理を実施する。

補正値算出部７は、例えば図２に示す補正値算出回路２５で実現されるものであり、輝度補正値算出部８、透過率補正値算出部９、画像再結合部１０及び補正値補間部１１を備えている。
補正値算出部７は、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する処理を実施する。
また、補正値算出部７は、大気透過率算出部６により波長帯毎に算出された各々の小領域における大気散乱の大気透過率から、パンクロマチック画像における分光放射輝度の補正値を算出する処理を実施する。
パンクロマチック画像は、マルチスペクトル画像が含んでいる複数の画像の波長帯を含む単一波長帯の画像である。

補正値算出部７の輝度補正値算出部８は、パラメータ記憶部２により記憶されている波長特性データを用いて、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、各々の小領域における分光放射輝度の第１の補正値をそれぞれ算出する処理を実施する。
補正値算出部７の透過率補正値算出部９は、パラメータ記憶部２により記憶されている波長特性データを用いて、大気透過率算出部６により波長帯毎に算出された各々の小領域における大気散乱の大気透過率から、各々の小領域における分光放射輝度の第２の補正値をそれぞれ算出する処理を実施する。

補正値算出部７の画像再結合部１０は、輝度補正値算出部８により算出された各々の小領域における分光放射輝度の第１の補正値を結合して、第１の補正値の２次元配列Ａを生成する処理を実施する。
また、画像再結合部１０は、透過率補正値算出部９により算出された各々の小領域における分光放射輝度の第２の補正値を結合して、第２の補正値の２次元配列Ｂを生成する処理を実施する。
２次元配列Ａ，Ｂにおける各々の要素は、画像分割部１により分割された各々の小領域に対応している。２次元配列Ａにおける各々の要素は、画像分割部１により分割された各々の小領域における分光放射輝度の第１の補正値を有しており、２次元配列Ｂにおける各々の要素は、画像分割部１により分割された各々の小領域における分光放射輝度の第２の補正値を有している。
補正値算出部７の補正値補間部１１は、画像再結合部１０により生成された２次元配列Ａ，Ｂの分解能が、パンクロマチック画像の分解能と一致するように、画像再結合部１０により生成された２次元配列Ａ，Ｂを補間する処理を実施する。

補正部１２は、図２に示す補正回路２６で実現されるものであり、放射輝度補正部１３及び透過率補正部１４を備えている。
補正部１２の放射輝度補正部１３は、パンクロマチックセンサからパンクロマチック画像を取得する処理を実施する。
放射輝度補正部１３は、補正値補間部１１により補間された２次元配列Ａにおける各々の要素が有する第１の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。
透過率補正部１４は、補正値補間部１１により補間された２次元配列Ｂにおける各々の要素が有する第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。

図１では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、補正値算出部７及び補正部１２のそれぞれが、図２に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、補正値算出回路２５及び補正回路２６で実現されるものを想定している。

ここで、パラメータ記憶回路２２は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）などの不揮発性又は揮発性の半導体メモリや、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｃ）などが該当する。
また、画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、補正値算出回路２５及び補正回路２６は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。

画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
ソフトウェア又はファームウェアはプログラムとして、コンピュータのメモリに格納される。コンピュータは、プログラムを実行するハードウェアを意味し、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）などが該当する。

図３は、画像処理装置がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。
画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２をコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、補正値算出部７及び補正部１２の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。
図４は、画像処理装置がソフトウェア又はファームウェアなどで実現される場合の処理手順を示すフローチャートである。

また、図２では、画像処理装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、画像処理装置がソフトウェアやファームウェアなどで実現される例を示しているが、画像処理装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。

次に動作について説明する。
この実施の形態１では、画像分割部１が、マルチスペクトルセンサから複数の波長帯の画像を含んでいるマルチスペクトル画像を取得し、放射輝度補正部１３が、パンクロマチックセンサからパンクロマチック画像を取得する。
このとき、マルチスペクトル画像が含んでいる複数の画像の波長帯のそれぞれは、帯域が狭く、パンクロマチック画像の波長帯は、マルチスペクトル画像が含んでいる複数の画像の波長帯を含んでいるものとする。
以下、マルチスペクトル画像が含んでいる複数の画像の波長帯をバンドｂと称し、この実施の形態１では、説明の便宜上、マルチスペクトル画像が、バンドｂ１の画像、バンドｂ２の画像、バンドｂ３の画像及びバンドｂ４の画像を含んでいるものとする。

画像分割部１は、マルチスペクトルセンサからマルチスペクトル画像を取得し（図４のステップＳＴ１）、図５に示すように、マルチスペクトル画像をＮ個の小領域Ｓ（ｉ，ｊ）に分割する（図４のステップＳＴ２）。
図５は、画像分割部１により分割されたマルチスペクトル画像を示す説明図である。
図５は、Ｎ＝１６の例を示しており、マルチスペクトル画像が１６個に分割されている。ｉ＝１，２，３，４、ｊ＝１，２，３，４である。
ただし、マルチスペクトル画像の分割数Ｎは、任意に設定できるものとし、例えば、マルチスペクトル画像に含まれている被写体の大きさなどから決定される。

具体的には、マルチスペクトル画像に含まれている被写体と、被写体の陰影などの暗部領域とが含まれる程度の小領域の大きさに分割されるように、マルチスペクトル画像の分割数Ｎを設定するようにする。
マルチスペクトル画像の中に、様々な大きさの被写体が含まれている場合には、マルチスペクトル画像に含まれている複数の被写体の中で、最も大きな被写体の大きさなどからマルチスペクトル画像の分割数Ｎを設定するようにして、被写体が含まれていない小領域Ｓ（ｉ，ｊ）の数を減らすようにしてもよい。

また、マルチスペクトル画像に靄又は霧などの濃度の変化を示す濃度分布が描画されている場合、図６に示すように、濃度分布において、濃度の極大値又は極小値を含む領域毎に、マルチスペクトル画像を分割するようにしてもよい。
図６は、画像分割部１により分割されたマルチスペクトル画像を示す説明図である。
図６では、図５と同様に、Ｎ＝１６の例を示しているが、Ｎ＜１６であって、濃度の極大値又は極小値を含む領域毎に分割された小領域Ｓ（ｉ，ｊ）が、図５のように分割された小領域Ｓ（ｉ，ｊ）よりも大きな領域であれば、図５のようにマルチスペクトル画像を分割する場合よりも、暗部輝度算出部４における暗部領域の分光放射輝度の算出数を低減することができる。
また、マルチスペクトル画像において、濃度の変化率が変化する箇所を把握できるように、マルチスペクトル画像が描画されている場合、小領域Ｓ（ｉ，ｊ）が、濃度の変化率が変化する箇所を含むようにマルチスペクトル画像を分割するようにしてもよい。

散乱輝度算出部３の暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）のヒストグラムをそれぞれ算出する（図４のステップＳＴ３）。
図７は、暗部輝度算出部４により算出されたバンドｂにおける小領域Ｓ（ｉ，ｊ）のヒストグラムを示す説明図である。
暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）の中で相対的に分光放射輝度が小さい暗部領域の分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）をそれぞれ算出する（図４のステップＳＴ４）。
図７では、小領域Ｓ（ｉ，ｊ）のヒストグラムの中で、分光放射輝度の最小値を、暗部領域の分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）としている例を示している。

ここでは、暗部輝度算出部４が、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出する例を示しているが、マルチスペクトル画像を構成している複数の画素の中に、画素値を有していない画素が含まれている場合がある。このような場合、暗部輝度算出部４は、マルチスペクトル画像を構成している複数の画素から、画素値を有していない画素を除いて、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出するようにしてもよい。
また、暗部輝度算出部４は、マルチスペクトル画像を構成している複数の画素の中に、画素値が常にゼロ付近の値になる欠陥画素が含まれている場合、マルチスペクトル画像を構成している複数の画素から、欠陥画素を除いて、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出するようにしてもよい。

散乱輝度算出部３の散乱輝度算出処理部５は、パラメータ記憶部２により記憶されている散乱特性データα（ｂ）を取得する（図４のステップＳＴ５）。
パラメータ記憶部２により記憶されている散乱特性データα（ｂ）は、バンドｂの散乱特性を示すデータであり、バンド固有の値である。
散乱特性データα（ｂ）は、例えば、ＭＯＤＴＲＡＮ（ＭＯＤｅｒａｔｅ ｒｅｓｏｌｕｔｉｏｎ ａｔｍｏｓｐｈｅｒｉｃ ＴＲＡＮｓｍｉｓｓｉｏｎ）などの光波大気伝搬計算シミュレータを用いて、各々の波長における大気散乱の放射輝度比から算出することができる。

散乱輝度算出処理部５は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を、以下の式（１）に示すように、散乱特性データα（ｂ）と、小領域Ｓ（ｉ，ｊ）の固有値Ｌ_{ｓｃａｔ０}（ｉ，ｊ）との積でモデル化する。
Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）＝α（ｂ）×Ｌ_{ｓｃａｔ０}（ｉ，ｊ） （１）
大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）は、バンドｂの大気散乱を示す分光放射輝度である。
小領域Ｓ（ｉ，ｊ）の固有値Ｌ_{ｓｃａｔ０}（ｉ，ｊ）は、波長依存性が無く、撮影環境に依存する小領域Ｓ（ｉ，ｊ）の固有値である。
図８は、散乱輝度算出処理部５によりモデル化された大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を示す説明図である。
図８において、実線は、大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）のモデルであり、〇は、小領域Ｓ（ｉ，ｊ）内の暗部領域の分光放射輝度である。

次に、散乱輝度算出処理部５は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、暗部領域の分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）と、式（１）の右辺であるα（ｂ）×Ｌ_{ｓｃａｔ０}（ｉ，ｊ）とが最も近くなるような固有値Ｌ_{ｓｃａｔ０}（ｉ，ｊ）を回帰分析によって推定する。
このとき、暗部領域の分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）の計測点は、マルチスペクトルセンサが有しているバンド数と同数である。また、回帰分析によるパラメータ決定手法の代表的な例として、最小二乗法などが挙げられる。

散乱輝度算出処理部５は、固有値Ｌ_{ｓｃａｔ０}（ｉ，ｊ）を推定すると、推定した固有値Ｌ_{ｓｃａｔ０}（ｉ，ｊ）と散乱特性データα（ｂ）を式（１）に代入して、バンドｂ毎に、各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）をそれぞれ算出する（図４のステップＳＴ６）。
散乱輝度算出処理部５は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を大気透過率算出部６及び輝度補正値算出部８に出力する。
ただし、散乱輝度算出処理部５は、ある小領域Ｓ（ｉ，ｊ）における暗部領域の分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）の計測点と、大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）のモデルとの相関が低い場合、当該小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を異常値として、例えば、負の値などを出力するようにする。
これにより、分光放射輝度Ｌ_Ｄ（ｂ，ｉ，ｊ）の計測点がモデルと類似していない暗部領域を含んでいる小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を除去することができる。
なお、計測点とモデルの相関が低い場合としては、例えば、最小二乗法における回帰分析の決定係数であるＲ二乗値が閾値Ｒｔｈ以下の場合などが挙げられる。

大気透過率算出部６は、パラメータ記憶部２により記憶されている散乱透過特性データｆ（ａ）を取得する（図４のステップＳＴ７）。ａ＝（ｂ，Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ））である。
パラメータ記憶部２により記憶されている散乱透過特性データｆ（ａ）は、小領域Ｓ（ｉ，ｊ）における大気散乱を示す大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）の大気透過率τ（ｂ，ｉ，ｊ）を示すバンドｂ毎の関数である。
散乱透過特性データｆ（ａ）は、例えば、ＭＯＤＴＲＡＮなどの光波大気伝搬計算シミュレータによって算出することができる。
図９は、バンドｂにおける小領域Ｓ（ｉ，ｊ）の大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）に対応する大気透過率τ（ｂ，ｉ，ｊ）を示す説明図である。

大気透過率算出部６は、以下の式（２）に示すように、散乱透過特性データｆ（ａ）の変数ａに、バンドｂと、散乱輝度算出処理部５から出力された小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）とを代入することで、バンドｂ毎に、小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を算出する（図４のステップＳＴ８）。
τ（ｂ，ｉ，ｊ）＝ｆ（ｂ，Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）） （２）
大気透過率算出部６は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を透過率補正値算出部９に出力する。
ただし、大気透過率算出部６は、ある小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）が想定値と数パーセント以上離れているような場合には、当該小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）を異常値として、例えば、負の値などを出力するようにする。

補正値算出部７の輝度補正値算出部８は、パラメータ記憶部２により記憶されている波長特性データを取得する（図４のステップＳＴ９）。
パラメータ記憶部２により記憶されている波長特性データは、マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯のそれぞれを示すデータである。
輝度補正値算出部８は、取得した波長特性データを用いて、散乱輝度算出処理部５から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）の重み係数ε（ｂ）を算出する。重み係数ε（ｂ）の算出処理の詳細は後述する。

輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）として、以下の式（３）に示すように、重み係数ε（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）の総和を算出する（図４のステップＳＴ１０）。

ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
Ｐａ’は、パンクロマチックセンサにより観測されるパンクロマチック画像に相当する画像の単一波長帯である。
輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を画像再結合部１０に出力する。
なお、輝度補正値算出部８は、散乱輝度算出処理部５から出力された小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）が異常値である場合、異常値である大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を用いずに、小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を算出するようにする。
例えば、輝度補正値算出部８は、大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）が異常値である小領域Ｓ（ｉ，ｊ）の周辺の小領域における大気散乱放射輝度を用いて、小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を補間する。
そして、輝度補正値算出部８は、補間した大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を用いて、第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を算出するようにする。

ここで、重み係数ε（ｂ）の算出処理について説明する。
図１０は、マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯を示す説明図である。
図１０の例では、パンクロマチックセンサの観測可能な波長帯は、λ_０（Ｐａ）〜λ_１（Ｐａ）である。
マルチスペクトルセンサの観測可能な波長帯は、λ_０（ｂ１）〜λ_１（ｂ１）、λ_０（ｂ２）〜λ_１（ｂ２）、λ_０（ｂ３）〜λ_１（ｂ３）及びλ_０（ｂ４）〜λ_１（ｂ４）である。

輝度補正値算出部８は、パンクロマチックセンサの観測可能な波長帯λ_０（Ｐａ）〜λ_１（Ｐａ）に対して、マルチスペクトルセンサの観測可能な波長帯λ’_０（ｂ）〜λ’_１（ｂ）を割り付けることで、以下の式（４）に示すように、重み係数ε（ｂ）を算出する。

ｂ＝ｂ１，ｂ２，ｂ３，ｂ４

ここで、輝度補正値算出部８におけるマルチスペクトルセンサの波長帯λ’_０（ｂ）〜λ’_１（ｂ）の割り付けは、例えば、波長帯λ’_０（ｂ）〜λ’_１（ｂ）が、マルチスペクトル画像が含んでいる画像のバンドｂ１，ｂ２，ｂ３，ｂ４のそれぞれと対応するように行う。
即ち、輝度補正値算出部８は、バンドｂ１に対して波長帯λ’_０（ｂ１）〜λ’_１（ｂ１）を割り付け、バンドｂ２に対して波長帯λ’_０（ｂ２）〜λ’_１（ｂ２）を割り付け、バンドｂ３に対して波長帯λ’_０（ｂ３）〜λ’_１（ｂ３）を割り付ける。
また、輝度補正値算出部８は、バンドｂ４に対して波長帯λ’_０（ｂ４）〜λ’_１（ｂ４）を割り付ける。
このような割り付けの場合、例えば、波長帯λ’_０（ｂ１），λ’_１（ｂ１）は、以下の式（５）及び式（６）のように表される。

重み係数ε（ｂ）における分子は、割り付けられた波長帯λ’_０（ｂ）〜λ’_１（ｂ）における分光放射輝度に対する波長積分を模擬するものである。
また、重み係数ε（ｂ）における分母は、パンクロマチックセンサの放射輝度から、分光放射輝度に変換すること意味している。
上記の式（３）は、割り付けられた波長帯λ’_０（ｂ）〜λ’_１（ｂ）における分光放射輝度を一定とみなしたものであり、パンクロマチックセンサの波長帯に対して、割り付けるバンドの数が多いほど、分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を高精度に算出することができる。

ここでは、輝度補正値算出部８が、式（４）によって、重み係数ε（ｂ）を算出する例を示しているが、以下の式（７）に示すように、パンクロマチックセンサにおける各波長λの分光感度Ｒ（λ）に基づいて重み係数ε（ｂ）を算出するようにしてもよい。この場合、パンクロマチックセンサが、各々のバンドｂの波長λで感度が著しく変化するようなセンサであっても、分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を高精度に算出することができる。

また、輝度補正値算出部８は、実画像のパンクロマチック画像を用いて、重み係数ε（ｂ）を算出するようにしてもよい。
例えば、輝度補正値算出部８は、実画像のパンクロマチック画像の分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）を取得する。
（ｘ，ｙ）は、パンクロマチック画像上の座標であり、マルチスペクトル画像における小領域の配列を示す（ｉ，ｊ）と異なる。
そして、輝度補正値算出部８は、式（３）から得られる分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を、マルチスペクトル画像から模擬したパンクロマチックセンサの分光放射輝度Ｌ（Ｐａ’，ｘ，ｙ）とみなし、式（３）において、分光放射輝度Ｌ（Ｐａ’，ｘ，ｙ）と分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）との差が最小になるような重み係数ε（ｂ）を算出する。

補正値算出部７の透過率補正値算出部９は、パラメータ記憶部２により記憶されている波長特性データを取得する（図４のステップＳＴ１１）。
透過率補正値算出部９は、取得した波長特性データを用いて、大気透過率算出部６から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）の重み係数ε’（ｂ）を算出する。重み係数ε’（ｂ）の算出処理の詳細は後述する。

透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）として、以下の式（８）に示すように、重み係数ε’（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）の総和を算出する（図４のステップＳＴ１２）。

ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を画像再結合部１０に出力する。
なお、透過率補正値算出部９は、大気透過率算出部６から出力された小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）が異常値である場合、異常値である大気透過率τ（ｂ，ｉ，ｊ）を用いずに、小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を算出するようにする。
例えば、透過率補正値算出部９は、大気透過率τ（ｂ，ｉ，ｊ）が異常値である小領域Ｓ（ｉ，ｊ）の周辺の小領域における大気透過率を用いて、小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を補間する。
そして、透過率補正値算出部９は、補間した大気透過率τ（ｂ，ｉ，ｊ）を用いて、第２の補正値τ（Ｐａ’，ｉ，ｊ）を算出するようにする。

ここで、重み係数ε’（ｂ）の算出処理について説明する。
大気透過率τ（ｂ，ｉ，ｊ）の重み係数ε’（ｂ）は、輝度補正値算出部８と同様に、式（４）によって算出することができる。
透過率補正値算出部９は、パンクロマチックセンサにおける各波長λの分光感度Ｒ（λ）及び光源の照度Ｅ（λ）を用いて、以下の式（９）に示すように、重み係数ε’（ｂ）を算出するようにしてもよい。

光源の照度Ｅ（λ）として、一般的な自然光である太陽の照度を用いることができる。
式（９）は、パンクロマチックセンサに入射される光源からの光のエネルギーに対する割り付け波長のエネルギー比を模擬している。

また、透過率補正値算出部９は、実画像のパンクロマチック画像を用いて、重み係数ε’（ｂ）を算出するようにしてもよい。
例えば、透過率補正値算出部９は、既知の分光放射輝度Ｌ_{ｓａｍｐｌｅ}を示す被写体と、暗部領域とが映っている実画像であるパンクロマチック画像を取得する。
そして、透過率補正値算出部９は、既知の分光放射輝度Ｌ_{ｓａｍｐｌｅ}と暗部領域の分光放射輝度との差分ｄＬを算出し、以下の式（１０）に示すように、差分ｄＬと既知の分光放射輝度Ｌ_{ｓａｍｐｌｅ}とから透過率τ（Ｐａ，ｘ，ｙ）を算出する。

透過率補正値算出部９は、式（８）から得られる第２の補正値τ（Ｐａ’，ｉ，ｊ）を、マルチスペクトル画像から模擬した透過率τ（Ｐａ’，ｘ，ｙ）とみなし、式（８）において、透過率τ（Ｐａ’，ｘ，ｙ）と透過率τ（Ｐａ，ｘ，ｙ）との差が最小になるような重み係数ε’（ｂ）を算出する。

補正値算出部７の画像再結合部１０は、輝度補正値算出部８から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を結合して、第１の補正値の２次元配列Ａを生成する（図４のステップＳＴ１３）。
また、画像再結合部１０は、透過率補正値算出部９から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を結合して、第２の補正値の２次元配列Ｂを生成する（図４のステップＳＴ１４）。
２次元配列Ａ，Ｂにおける各々の要素（ｉ，ｊ）は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。マルチスペクトル画像の分割数がＮで、小領域Ｓ（ｉ，ｊ）の個数がＮであれば、２次元配列の要素（ｉ，ｊ）の数もＮである。
２次元配列Ａにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を有しており、２次元配列Ｂにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を有している。
画像再結合部１０は、生成した２次元配列Ａ，Ｂを補正値補間部１１に出力する。

補正値算出部７の補正値補間部１１は、画像再結合部１０から出力された２次元配列Ａ，Ｂの分解能のそれぞれが、パンクロマチック画像の分解能と一致するように、画像再結合部１０から出力された２次元配列Ａ，Ｂのそれぞれを補間する（図４のステップＳＴ１５）。
例えば、パンクロマチック画像のｘ方向の画素数が、２次元配列Ａ，Ｂにおけるｘ方向の要素（ｉ，ｊ）の数のＭｘ倍であれば、２次元配列Ａ，Ｂのｘ方向をＭｘ倍にアップサンプリングする。
また、パンクロマチック画像のｙ方向の画素数が、２次元配列Ａ，Ｂにおけるｙ方向の要素（ｉ，ｊ）の数のＭｙ倍である場合、２次元配列Ａ，Ｂのｙ方向をＭｙ倍にアップサンプリングする。
これにより、２次元配列Ａ，Ｂの分解能が、パンクロマチック画像の分解能と同等になる。

２次元配列Ａ，Ｂの補間処理としては、２次元配列Ａ，Ｂの要素間の補正値が連続的となるような補間処理が望ましく、例えば、バイリニア又はＣｕｂｉｃ Ｃｏｎｖｏｌｕｔｉｏｎなどの補間処理を実施することが考えられる。また、２次元配列Ａ，Ｂの要素間の補正値が連続的かつ曲線を描くように、ベジエ曲線又はＳｐｌｉｎｅ曲線などの補間処理を実施するようにしてもよい。
これにより、第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）は、第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｘ，ｙ）が補正値補間部１１から放射輝度補正部１３に出力される。また、第２の補正値τ（Ｐａ’，ｉ，ｊ）は、第２の補正値τ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第２の補正値τ（Ｐａ’，ｘ，ｙ）が補正値補間部１１から透過率補正部１４に出力される。
なお、２次元配列Ａの或る要素が有している第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）が異常値である場合、補正値補間部１１は、当該要素の周囲に存在している要素が有している補正値に基づいて、当該要素が有している補正値を補間してから、２次元配列Ａの分解能を変換する。
また、２次元配列Ｂの或る要素が有している第２の補正値τ（Ｐａ’，ｉ，ｊ）が異常値である場合、補正値補間部１１は、当該要素の周囲に存在している要素が有している補正値に基づいて、当該要素が有している補正値を補間してから、２次元配列Ｂの分解能を変換する。

補正部１２の放射輝度補正部１３は、パンクロマチックセンサからパンクロマチック画像を取得する（図４のステップＳＴ１６）。
放射輝度補正部１３は、以下の式（１１）に示すように、取得したパンクロマチック画像の分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）から、補正値補間部１１から出力された第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｘ，ｙ）を減算することで、パンクロマチック画像における分光放射輝度を補正する（図４のステップＳＴ１７）。
Ｌ１（ｘ，ｙ）
＝Ｌ（Ｐａ，ｘ，ｙ）−Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ） （１１）
式（１１）において、Ｌ１（ｘ，ｙ）は、補正後の分光放射輝度である。
放射輝度補正部１３は、分光放射輝度を補正したパンクロマチック画像を透過率補正部１４に出力する。

補正部１２の透過率補正部１４は、以下の式（１２）に示すように、放射輝度補正部１３から出力されたパンクロマチック画像の補正後の分光放射輝度Ｌ１（ｘ，ｙ）を、補正値補間部１１から出力された第２の補正値τ（Ｐａ’，ｘ，ｙ）で除算することで、さらに、パンクロマチック画像における分光放射輝度を補正する（図４のステップＳＴ１８）。

式（１２）において、Ｌ２（ｘ，ｙ）は、補正後の分光放射輝度である。
透過率補正部１４は、分光放射輝度を補正したパンクロマチック画像を出力する。

以上で明らかなように、この実施の形態１によれば、散乱輝度算出部３により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部７を設け、補正部１２が、補正値算出部７により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、パンクロマチックセンサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる効果を奏する。

また、この実施の形態１によれば、補正値算出部７が、大気透過率算出部６により波長帯毎に算出された各々の領域における大気透過率から、パンクロマチック画像における分光放射輝度の第２の補正値を算出し、補正部１２が、第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、第１の補正値だけを用いて、パンクロマチック画像における分光放射輝度を補正する場合よりも、コントラストの補正精度を高めることができる効果を奏する。

実施の形態２．
上記実施の形態１では、補正部１２が、第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｘ，ｙ）及び第２の補正値τ（Ｐａ’，ｘ，ｙ）を用いて、パンクロマチック画像における分光放射輝度を補正する例を示している。
この実施の形態２では、パンクロマチック画像に含まれているオフセット成分を除去するために第３の補正値を算出し、第３の補正値を用いて、さらに、パンクロマチック画像における分光放射輝度を補正する例を説明する。

図１１は、この発明の実施の形態２による画像処理装置を示す構成図であり、図１２は、この発明の実施の形態２による画像処理装置を示すハードウェア構成図である。
図１１及び図１２において、図１及び図２と同一符号は同一または相当部分を示すので説明を省略する。
反射輝度算出部４１は、例えば図１２に示す反射輝度算出回路５１で実現される。
反射輝度算出部４１は、マルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）と、散乱輝度算出処理部５によりバンドｂ毎に算出された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）とから、各々の小領域Ｓ（ｉ，ｊ）における大気反射を示す大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）をそれぞれ算出する処理を実施する。

補正値算出部４２は、例えば図２に示す補正値算出回路５２で実現されるものであり、輝度補正値算出部８、透過率補正値算出部９、反射補正値算出部４３、画像再結合部４４及び補正値補間部４５を備えている。
補正値算出部４２の反射補正値算出部４３は、反射輝度算出部４１によりバンドｂ毎に算出された各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）から、各々の小領域における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）をそれぞれ算出する処理を実施する。

補正値算出部４２の画像再結合部４４は、上記実施の形態１の画像再結合部１０と同様に、２次元配列Ａ及び２次元配列Ｂを生成する処理を実施する。
また、画像再結合部４４は、反射補正値算出部４３により算出された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を結合して、第３の補正値の２次元配列Ｃを生成する処理を実施する。
２次元配列Ｃにおける各々の要素は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。２次元配列Ｃにおける各々の要素は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値を有している。
補正値算出部４２の補正値補間部４５は、画像再結合部４４により生成された２次元配列Ａ，Ｂ，Ｃの分解能が、パンクロマチック画像の分解能と一致するように、画像再結合部４４により生成された２次元配列Ａ，Ｂ，Ｃを補間する処理を実施する。

補正部４６は、図１２に示す補正回路５３で実現されるものであり、放射輝度補正部１３、透過率補正部１４及び反射補正部４７を備えている。
補正部４６の反射補正部４７は、補正値補間部４５により補間された２次元配列Ｃにおける各々の要素が有する第３の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。

図１１では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２及び補正部４６のそれぞれが、図１２に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２及び補正回路５３で実現されるものを想定している。

画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２及び補正回路５３は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、または、これらを組み合わせたものが該当する。
画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。

画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２を図３に示すコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２及び補正部４６の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。
また、図１２では、画像処理装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、画像処理装置がソフトウェアやファームウェアなどで実現される例を示しているが、画像処理装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。

次に動作について説明する。
靄又は霧などの濃度が高い場合、被写体とセンサの間に存在している大気中の粒子に光が散乱される大気散乱のほかに、大気中の粒子に光が反射される大気反射が生じることがある。
大気反射が生じることで、大気からの反射光成分がオフセット成分としてマルチスペクトル画像及びパンクロマチック画像に重畳される。
この実施の形態２では、補正部４６の反射補正部４７が、パンクロマチック画像に重畳されているオフセット成分を除去する処理を実施する。
以下、上記実施の形態１と相違する箇所を説明する。

反射輝度算出部４１は、画像分割部１により分割されたマルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）を取得する。
反射輝度算出部４１は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）を、以下の式（１３）に示すように、散乱輝度算出処理部５により算出された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）と、各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）との和でモデル化する。
Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）
＝Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）＋Ｌｒ（ｂ，ｉ，ｊ） （１３）
図１３は、反射輝度算出部４１によりモデル化された分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）を示す説明図である。
図１３において、実線は、分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）のモデルであり、〇は、小領域Ｓ（ｉ，ｊ）内の暗部領域の分光放射輝度である。

ここで、大気からの反射光成分の分光特性が変化せずに、反射光成分の絶対値のみが変化するとすれば、マルチスペクトル画像の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）は、以下の式（１４）のように表すことができる。
Ｌｒ（ｂ，ｉ，ｊ）＝β（ｂ）Ｌ_ｒ０（ｉ，ｊ） （１４）
式（１４）において、β（ｂ）は、反射光の分光特性成分、Ｌ_ｒ０（ｉ，ｊ）は、反射光成分の固有値である。
反射光の分光特性成分β（ｂ）は、以下の式（１５）に示すように、白色の被写体が持つ分光放射輝度より算出することができる。

式（１５）において、ρ（λ）は、白色の被写体の反射率、Ｅ（λ）は、光源である太陽の照度である。
式（１３）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）は、式（１）のように表されるので、式（１３）は、以下の式（１６）のように表される。
Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）
＝α（ｂ）Ｌ_{ｓｃａｔ０}（ｉ，ｊ）＋β（ｂ）Ｌ_ｒ０（ｉ，ｊ） （１６）

反射輝度算出部４１は、式（１６）のモデルにおいて、反射光成分の固有値Ｌ_ｒ０（ｉ，ｊ）を回帰分析によって推定する。
反射輝度算出部４１は、反射光成分の固有値Ｌ_ｒ０（ｉ，ｊ）を推定すると、推定した固有値Ｌ_ｒ０（ｉ，ｊ）と、反射光の分光特性成分β（ｂ）とを式（１４）に代入して、マルチスペクトル画像の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を算出する。
反射輝度算出部４１は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を反射補正値算出部４３に出力する。

補正値算出部４２の反射補正値算出部４３は、パラメータ記憶部２により記憶されている波長特性データを取得する。
反射補正値算出部４３は、取得した波長特性データを用いて、反射輝度算出部４１から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）の重み係数ε”（ｂ）を算出する。重み係数ε”（ｂ）は、例えば、重み係数ε（ｂ）又は重み係数ε’（ｂ）と同様に算出することができる。

反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）として、以下の式（１７）に示すように、重み係数ε”（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）の総和を算出する。

ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力する。

補正値算出部４２の画像再結合部４４は、上記実施の形態１の画像再結合部１０と同様に、第１の補正値の２次元配列Ａ及び第２の補正値の２次元配列Ｂのそれぞれを生成する。
また、画像再結合部４４は、反射補正値算出部４３から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を結合して、第３の補正値の２次元配列Ｃを生成する。
２次元配列Ｃにおける各々の要素（ｉ，ｊ）は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。マルチスペクトル画像の分割数がＮで、小領域Ｓ（ｉ，ｊ）の個数がＮであれば、２次元配列Ｃの要素（ｉ，ｊ）の数もＮである。
２次元配列Ｃにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を有している。
画像再結合部４４は、生成した２次元配列Ａ，Ｂ，Ｃを補正値補間部４５に出力する。

補正値算出部４２の補正値補間部４５は、上記実施の形態１の補正値補間部１１と同様に、画像再結合部４４から出力された２次元配列Ａ，Ｂのそれぞれを補間する。
また、補正値補間部４５は、画像再結合部４４から出力された２次元配列Ｃを補間する。
例えば、パンクロマチック画像のｘ方向の画素数が、２次元配列Ｃにおけるｘ方向の要素（ｉ，ｊ）の数のＭｘ倍であれば、２次元配列Ｃのｘ方向をＭｘ倍にアップサンプリングする。
これにより、第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）は、第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）が補正値補間部４５から反射補正部４７に出力される。

補正部４６の反射補正部４７は、以下の式（１８）に示すように、放射輝度補正部１３から出力されたパンクロマチック画像の補正後の分光放射輝度Ｌ１（ｘ，ｙ）から、補正値補間部４５から出力された第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）を減算することで、パンクロマチック画像における分光放射輝度を補正する。
Ｌ３（ｘ，ｙ）＝Ｌ１（ｘ，ｙ）−Ｌｒ（Ｐａ’，ｘ，ｙ） （１８）
式（１８）において、Ｌ３（ｘ，ｙ）は、補正後の分光放射輝度である。
反射補正部４７は、分光放射輝度を補正したパンクロマチック画像を透過率補正部１４に出力する。
この実施の形態２では、反射補正部４７が放射輝度補正部１３の後段に設けられている例を示しているが、反射補正部４７が放射輝度補正部１３の前段に設けられていてもよい。

以上で明らかなように、この実施の形態２によれば、補正値算出部４２が、反射輝度算出部４１により波長帯毎に算出された各々の領域における大気反射放射輝度から、パンクロマチック画像における分光放射輝度の第３の補正値を算出し、補正部４６が、補正値算出部４２により算出された第１及び第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するほかに、補正値算出部４２により算出された第３の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、上記実施の形態１，２と同様の効果を奏するほかに、マルチスペクトル画像及びパンクロマチック画像に重畳されている大気からの反射光成分を除去することができる効果を奏する。

実施の形態３．
この実施の形態３では、マルチスペクトル画像における分光放射輝度についても補正し、分光放射輝度が補正されたパンクロマチック画像と、分光放射輝度が補正されたマルチスペクトル画像との色合成を行う例を説明する。

図１４は、この発明の実施の形態３による画像処理装置を示す構成図であり、図１５は、この発明の実施の形態３による画像処理装置を示すハードウェア構成図である。
図１４及び図１５において、図１１及び図１２と同一符号は同一または相当部分を示すので説明を省略する。
この実施の形態３では、補正値算出部４２は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度の補正値を算出するほかに、マルチスペクトル画像における分光放射輝度の補正値を算出する処理を実施する。
この実施の形態３では、補正部４６は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度を補正するほかに、マルチスペクトル画像における分光放射輝度を補正する処理を実施する。

色合成処理部４８は、図１５に示す色合成処理回路５４で実現される。
色合成処理部４８は、補正部４６により分光放射輝度が補正されたパンクロマチック画像と補正部４６により分光放射輝度が補正されたマルチスペクトル画像との色合成処理を実施する。
図１４では、色合成処理部４８が、上記実施の形態２における図１１の画像処理装置に適用されている例を示しているが、色合成処理部４８が、上記実施の形態１における図１の画像処理装置に適用されているものであってもよい。

図１４では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２、補正部４６及び色合成処理部４８のそれぞれが、図１５に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２、補正回路５３及び色合成処理回路５４で実現されるものを想定している。

画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２、補正回路５３及び色合成処理回路５４は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、または、これらを組み合わせたものが該当する。
画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。

画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２を図３に示すコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２、補正部４６及び色合成処理部４８の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。

次に動作について説明する。
以下、上記実施の形態１，２と相違する箇所を説明する。
この実施の形態３では、輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌ_ｓｃａｔ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌ_ｓｃａｔ（ｂ，ｉ，ｊ）を第４の補正値として画像再結合部４４に出力する。
この実施の形態３では、透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）を第５の補正値として画像再結合部４４に出力する。
この実施の形態３では、反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）として画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を第６の補正値として画像再結合部４４に出力する。

補正値算出部４２の画像再結合部４４は、上記実施の形態２と同様に、第１の補正値の２次元配列Ａを生成するほかに、輝度補正値算出部８から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第４の補正値であるＬ_ｓｃａｔ（ｂ，ｉ，ｊ）を結合して、第４の補正値の２次元配列Ｄを生成する。
また、画像再結合部４４は、第２の補正値の２次元配列Ｂを生成するほかに、透過率補正値算出部９から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第５の補正値である大気透過率τ（ｂ，ｉ，ｊ）を結合して、第５の補正値の２次元配列Ｅを生成する。
また、画像再結合部４４は、第３の補正値の２次元配列Ｃを生成するほかに、反射補正値算出部４３から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第６の補正値である大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を結合して、第６の補正値の２次元配列Ｆを生成する。

補正値算出部４２の補正値補間部４５は、上記実施の形態２と同様に、画像再結合部４４から出力された２次元配列Ａ，Ｂ，Ｃのそれぞれを補間する。
また、補正値補間部４５は、画像再結合部４４から出力された２次元配列Ｄ，Ｅ，Ｆのそれぞれを補間する。
これにより、第４の補正値であるＬ_ｓｃａｔ（ｂ，ｉ，ｊ）は、Ｌ_ｓｃａｔ（ｂ，ｘ，ｙ）に変換され、変換後の第４の補正値であるＬ_ｓｃａｔ（ｂ，ｘ，ｙ）が補正値補間部４５から放射輝度補正部１３に出力される。
第５の補正値である大気透過率τ（ｂ，ｉ，ｊ）は、τ（ｂ，ｘ，ｙ）に変換され、変換後の第５の補正値であるτ（ｂ，ｘ，ｙ）が補正値補間部４５から透過率補正部１４に出力される。
第６の補正値である大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）は、Ｌｒ（ｂ，ｘ，ｙ）に変換され、変換後の第６の補正値であるＬｒ（ｂ，ｘ，ｙ）が補正値補間部４５から反射補正部４７に出力される。

放射輝度補正部１３は、上記実施の形態１，２と同様に、パンクロマチック画像における分光放射輝度を補正する。
また、放射輝度補正部１３は、マルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）から、補正値補間部４５から出力された第４の補正値であるＬ_ｓｃａｔ（ｂ，ｘ，ｙ）を減算することで、マルチスペクトル画像における分光放射輝度を補正する。

反射補正部４７は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度を補正する。
また、反射補正部４７は、放射輝度補正部１３から出力されたマルチスペクトル画像の補正後の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）から、補正値補間部４５から出力された第６の補正値であるＬｒ（ｂ，ｘ，ｙ）を減算することで、マルチスペクトル画像における分光放射輝度を補正する。

透過率補正部１４は、上記実施の形態１，２と同様に、パンクロマチック画像における分光放射輝度を補正する。
また、透過率補正部１４は、反射補正部４７から出力されたマルチスペクトル画像の補正後の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌ_{ｓｅｎｓｏｒ}（ｂ，ｉ，ｊ）を、補正値補間部４５から出力された第５の補正値であるτ（ｂ，ｘ，ｙ）で除算することで、マルチスペクトル画像における分光放射輝度を補正する。

色合成処理部４８は、補正部４６により分光放射輝度が補正されたパンクロマチック画像と、補正部４６により分光放射輝度が補正されたマルチスペクトル画像との色合成処理を実施し、色合成処理後の画像を出力する。
色合成処理として、例えば、ＩＨＳ（Ｉｎｔｅｎｓｉｔｙ Ｈｕｅ Ｓａｔｕｒａｔｉｏｎ）変換と呼ばれる色合成処理を用いることができる。ＩＨＳ変換と呼ばれる色合成処理では、マルチスペクトル画像から色空間である輝度、色相及び彩度を算出し、この算出した輝度をパンクロマチック画像に置き換えることで、パンクロマチック画像の分解能を持つ輝度情報に、色相及び彩度の情報を加えるものである。
色合成処理としては、ｂｒｏｖｅｙ変換、Ｇｒａｍ−Ｓｃｈｍｉｄｔ変換、あるいは、主成分分析を用いた変換処理なども使用することができる。
高解像度のパンクロマチック画像と低解像度のマルチスペクトル画像とを合成し、高解像度かつ色情報を保有する画像は、パンシャープン画像と呼ばれることがある。

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

この発明は、パンクロマチック画像における分光放射輝度を補正する画像処理装置に適している。

１ 画像分割部、２ パラメータ記憶部、３ 散乱輝度算出部、４ 暗部輝度算出部、５ 散乱輝度算出処理部、６ 大気透過率算出部、７ 補正値算出部、８ 輝度補正値算出部、９ 透過率補正値算出部、１０ 画像再結合部、１１ 補正値補間部、１２ 補正部、１３ 放射輝度補正部、１４ 透過率補正部、２１ 画像分割回路、２２ パラメータ記憶回路、２３ 散乱輝度算出回路、２４ 透過率算出回路、２５ 補正値算出回路、２６ 補正回路、３１ メモリ、３２ プロセッサ、４１ 反射輝度算出部、４２ 補正値算出部、４３ 反射補正値算出部、４４ 画像再結合部、４５ 補正値補間部、４６ 補正部、４７ 反射補正部、４８ 色合成処理部、５１ 反射輝度算出回路、５２ 補正値算出回路、５３ 補正回路、５４ 色合成処理回路。

Claims (4)

1. 複数の波長帯の画像を含んでいるマルチスペクトル画像を複数の領域に分割する画像分割部と、
   前記マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、前記画像分割部により分割された各々の領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する散乱輝度算出部と、
   前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部と、
   前記補正値算出部により算出された分光放射輝度の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正する補正部と
   を備えた画像処理装置。
2. 前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記波長帯毎に、前記画像分割部により分割された各々の領域における大気散乱の大気透過率をそれぞれ算出する大気透過率算出部を備え、
   前記補正値算出部は、前記パンクロマチック画像における分光放射輝度の補正値として、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記パンクロマチック画像における分光放射輝度の第１の補正値を算出するほかに、前記大気透過率算出部により波長帯毎に算出された各々の領域における大気散乱の大気透過率から、前記パンクロマチック画像における分光放射輝度の第２の補正値を算出し、
   前記補正部は、前記補正値算出部により算出された第１の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正するほかに、前記補正値算出部により算出された第２の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正することを特徴とする請求項１記載の画像処理装置。
3. 前記マルチスペクトル画像の波長帯毎の各々の領域における分光放射輝度と、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度とから、前記波長帯毎に、前記画像分割部により分割された各々の領域における大気反射を示す大気反射放射輝度をそれぞれ算出する反射輝度算出部を備え、
   前記補正値算出部は、前記第１及び第２の補正値を算出するほかに、前記反射輝度算出部により波長帯毎に算出された各々の領域における大気反射放射輝度から、前記パンクロマチック画像における分光放射輝度の第３の補正値を算出し、
   前記補正部は、前記補正値算出部により算出された第１及び第２の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正するほかに、前記補正値算出部により算出された第３の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正することを特徴とする請求項２記載の画像処理装置。
4. 前記補正値算出部は、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記波長帯毎に、前記マルチスペクトル画像における分光放射輝度の補正値を算出し、
   前記補正部は、前記補正値算出部により算出されたマルチスペクトル画像における分光放射輝度の補正値を用いて、前記マルチスペクトル画像における分光放射輝度を補正し、
   前記補正部により分光放射輝度が補正されたパンクロマチック画像と前記補正部により分光放射輝度が補正されたマルチスペクトル画像との色合成を行う色合成処理部を備えたことを特徴とする請求項１記載の画像処理装置。

Images