(57)【要約】【課題】ミクセルを含む画像データから、周期的な空間分布を有する対象物の占有率を推定する。【解決手段】画像解析装置は、周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号を含む画像データから、対象物の占有率を算出するものであり、分析部と、基準値設定部と、任意性消去部とを有する。分析部は、対象物の反射率と背景の反射率との差を混合比とすると共に、対象物の占有率を独立成分として、独立成分分析により混合信号からノイズを分離して、混合比及び独立成分を任意性のある解として算出する。基準値設定部は、背景の反射率の波長微分において安定した値を、基準値として求めておく。任意性消去部は、背景の反射率が基準値をとると仮定して、任意性のある解として算出された混合比及び独立成分を、それぞれ１つの解に決定する。従って、前述の混合信号から、対象物の占有率を算出できる。【選択図】      図１
2004184391.tif 000002
【特許請求の範囲】
【請求項１】
周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号を含む画像データを解析する画像解析装置であって、
独立成分分析における混合比を『前記対象物の反射率』と『前記背景の反射率』との差とすると共に、独立成分分析における独立成分を『前記対象物の占有率』として、独立成分分析により前記混合信号からノイズを分離して、前記混合比及び前記独立成分を任意性のある解として算出する分析部と、
前記背景の反射率の波長微分値において、前記背景の種類に拘わらずに同じ値とみなせるものを、基準値として求めておく基準値設定部と、
前記背景の反射率が前記基準値をとると仮定して、任意性のある解として算出された前記混合比及び前記独立成分を、それぞれ１つの解に決定する任意性消去部とを備えていることを特徴とする画像解析装置。
【請求項２】
請求項１記載の画像解析装置において、
前記分析部は、
前記画像データが撮像されたときの入射光強度をI0(λ)、ノイズ成分をN(λ，x)、対象物の占有率をηV(x)、背景の占有率をηS(x)、対象物の反射率をRV(λ)、背景の反射率をRS(λ)として（ここで、(x)は前記画像データにおける画素の位置ｘの関数であることを示し、(λ)は波長λの関数であることを示し、(λ，x)はｘ及びλの関数であることを示す。）、
ηS(x)＝１−ηV(x)  とみなし、次式
2004184391.tif 000014
の関係に基づいて、RV(λ)とRS(λ)との差を前記混合比として独立成分分析を行うことを特徴とする画像解析装置。
【請求項３】
請求項１または請求項２記載の画像解析装置において、
前記基準値設定部は、
前記対象物の反射率の微分値がゼロとなる波長λi(i＝１,２,・・・,ｐ、ｐは２以上の自然数) における、前記背景の反射率値からなる行列をＲｅｆSJとして、
前記ＲｅｆSJの波長微分で与えられる行列ＲｅｆSJ’と、前記ＲｅｆSJ’のミラー行列とを判別分析により写像して、写像後の行列の要素から前記基準値を選択することを特徴とする画像解析装置。
【請求項４】
請求項１〜請求項３のいずれか１項記載の前記分析部、前記基準値設定部、及び前記任意性消去部として、コンピュータを機能させるための画像解析プログラム。
【請求項５】
周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号を含む画像データを解析する画像解析方法であって、
独立成分分析における混合比を『前記対象物の反射率』と『前記背景の反射率』との差とすると共に、独立成分分析における独立成分を『前記対象物の占有率』として、独立成分分析により前記混合信号からノイズを分離して、前記混合比及び前記独立成分を任意性のある解として算出する手順と、
前記背景の反射率の波長微分値において、前記背景の種類に拘わらずに同じ値とみなせるものを、基準値として求めておく手順と、
前記背景の反射率が前記基準値をとると仮定して、任意性のある解として算出された前記混合比及び前記独立成分を、それぞれ１つの解に決定する手順とを有することを特徴とする画像解析方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、画像解析装置、画像解析プログラム、及び画像解析方法に関する。具体的には、上空から撮影された土地の画像データなどが有する混合スペクトルを複数のピュアスペクトルに分離して、その被覆パターンを解析する技術に関する。
【０００２】
【従来の技術】
リモートセンシングは、航空機画像等を用いることで、立ち入り困難な地域の地理的状況を、調査に赴く事なくある程度把握できるという利点を有している。このため、リモートセンシングは、農業、環境、防災等の様々な分野においての利用が期待されている。
【０００３】
通常、上空から撮影した衛星画像や航空機映像では、撮像系の解像度に限界があるため、対象とする構造物と背景とが、１画素内に複数混在することが多い。即ち、画像データにおける少なくとも一部の画素値（スペクトル）は、その画素内における各対象物の値（ピュアスペクトル）を、混合した値（混合スペクトル）となっている。このように混合スペクトルとなっている画素を、ミクセルという。
【０００４】
ミクセルを含む画像データから被覆パターン（各対象物の占有率）を求めるためには、混合スペクトルを複数のピュアスペクトルに分離する作業を、全てのミクセルに対して行う必要がある。
ミクセルを含む画像データを解析する場合、従来の手法では、混合スペクトルは既知な構成要素のピュアスペクトルの線形和であると仮定している（例えば、非特許文献１参照）。しかし、このように構成要素を全て既知として扱うことは現実的でなく、各構成要素もそれらの混合比も未知である場合が殆どである。そして、各構成要素が推測不可能な場合に混合スペクトルをピュアスペクトルに分離することは、極めて困難であるため、極く一部の分野でしか利用可能な精度で実現できていない。
【０００５】
混合スペクトルをピュアスペクトルに分離する方法が既に実用化されている分野の１つとして、ハイパースペクトル計測を用いた鉱物等の資源探査の分野が知られている（例えば、非特許文献２、及び非特許文献３参照）。なお、ハイパースペクトル計測は、数ナノメートルの波長分解能で、対象物の反射スペクトルを計測するものであり、これにより得られる観測データをハイパースペクトルデータという。また、ここでの反射スペクトルは、数十〜数百種類の波長に対して輝度値を有するハイパースペクトルとなる。
【０００６】
上述した鉱物の資源探査方法は、ハイパースペクトル計測の高波長分解能という特徴と、鉱物毎に局所的な吸収特性が異なることとを利用したものである。具体的には、線形モデルを多次元空間に展開し、スペクトル空間における局所的なピークをピュアスペクトルとして抽出する。この後、鉱物毎の吸収特性との整合性をとることにより、高精度に特定の鉱物を検出する。
【０００７】
【非特許文献１】
著者：小熊宏之、島崎彦人、山形与志樹著
題名：『森林樹冠率推定におけるミクセル分類手法の有効性に関する研究』
出典：写真測量とリモートセンシング、Vol.41、No2、p.4-14
出版社：写真測量学会  2002年4月出版
【非特許文献２】
著者：J.W.Boardman
題名：『AUTOMATING SPECTRAL UNMIXING OF AVIRIS DATA USING CONVEXGEOMETRY CONCEPTS』
出典：Fourth JPL Airborne Geoscience Workshop, JPLPublication 93-26, Vol.1, pp. 11-14 (1993)
http://makalu.jpl.nasa.gov/docs/workshops/93#docs/4.PDF
【非特許文献３】
著者：J.W.Boardman, F.A.Kruse and R.O.Green
題名：『MAPPING TARGET SIGNATURES VIA PARTIAL UNMIXINGOF AVIRIS DATA』
出典：Fifth JPL Airborne Geoscience WorkshopJPL Publication 95-1, Vol.1, pp. 23-26 (1995)
http://makalu.jpl.nasa.gov/docs/workshops/95#docs/7.PDF
【０００８】
【発明が解決しようとする課題】
上述のように、異なる局所的な吸収特性を有する物質のピュアスペクトルは、混合スペクトルから分離可能である。しかし、局所的な吸収特性を有さない物質や、マルチスペクトル（数種類〜数十種類の波長に対して輝度値を有するスペクトル）のような低波長分解能のスペクトルには、同様の手法を適用できない。
【０００９】
本発明の目的は、ミクセルを含む画像データから、周期的な空間分布を有する対象物の占有率分布を高精度に推定する技術を提供することである。
【００１０】
【課題を解決するための手段】
本発明者は、各構成要素が未知である混合スペクトルからピュアスペクトルを抽出する手法として、独立成分分析による周期的ハイパースペクトルデータの分離方法を適用することを考えた。ここで、独立成分分析（Independent Component Analysis、以下、ＩＣＡと略記）は、複数の音源が発する音を複数の場所で観測して、複数の観測信号に基づいて、個々の音源が発した原信号をそれぞれ推定する技術である。具体的には、
複数の観測信号に対応する混合信号ベクトルをz＝(z1,・・・,zK)T 、
原信号ベクトルをs＝(s1,・・・,sU)T、
混合比が要素aij を持つK行U列の行列をA 、
無相関のガウス性ノイズ成分をz0
とする。このとき、次式の関係が成り立つ。
【数１】
2004184391.tif 000003
ＩＣＡは、観測したzからガウス性ノイズ成分z0を分離して、混合比Aと原信号ベクトルsとを同時に推定する手法である。このＩＣＡを利用して、本発明は、以下のように構成される。
【００１１】
＜請求項１＞
請求項１の画像解析装置は、周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号を含む画像データを解析するものであり、分析部と、基準値設定部と、任意性消去部とを備えていることを特徴とする。分析部は、独立成分分析における混合比を対象物の反射率と背景の反射率との差とすると共に、独立成分分析における独立成分を対象物の占有率として、独立成分分析により混合信号からノイズを分離して、混合比及び独立成分を任意性のある解として算出する。基準値設定部は、背景の反射率の波長微分値において、背景の種類に拘わらずに同じ値とみなせるものを、基準値として求めておく。任意性消去部は、背景の反射率が基準値をとると仮定して、任意性のある解として算出された混合比及び独立成分を、それぞれ１つの解に決定する。
【００１２】
＜請求項２＞
請求項２の画像解析装置は、請求項１記載の発明において、以下の点を特徴とする。第１に、分析部は、画像データが撮像されたときの入射光強度をI0(λ)、ノイズ成分をN(λ，x)、対象物の占有率をηV(x)、背景の占有率をηS(x)、対象物の反射率をRV(λ)、背景の反射率をRS(λ)として（ここで、(x)は画像データにおける画素の位置ｘの関数であることを示し、(λ)は波長λの関数であることを示し、(λ，x)はｘ及びλの関数であることを示す。）、
ηS(x)＝１−ηV(x)  とみなす。
第２に、分析部は、次式
2004184391.tif 000004
の関係に基づいて、RV(λ)とRS(λ)との差を混合比として独立成分分析を行う。
【００１３】
＜請求項３＞
請求項３の画像解析装置は、請求項１または請求項２記載の発明において、『基準値設定部は、対象物の反射率の微分値がゼロとなる波長λi(i＝１,２,・・・,ｐ、ｐは２以上の自然数) における、背景の反射率値からなる行列をＲｅｆSJとして、ＲｅｆSJの波長微分で与えられる行列ＲｅｆSJ’と、ＲｅｆSJ’のミラー行列とを判別分析により写像して、写像後の行列の要素から基準値を選択する』ことを特徴とする。
【００１４】
＜請求項４＞
請求項４の画像解析プログラムは、請求項１〜請求項３のいずれか１項記載の分析部、基準値設定部、及び任意性消去部として、コンピュータを機能させるための画像解析プログラムである。
【００１５】
＜請求項５＞
請求項５の画像解析方法は、周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号を含む画像データを解析するものであって、以下の手順を有することを特徴とする。
１つは、独立成分分析における混合比を『対象物の反射率』と『背景の反射率』との差とすると共に、独立成分分析における独立成分を『対象物の占有率』として、独立成分分析により混合信号からノイズを分離して、混合比及び独立成分を任意性のある解として算出する手順である。
１つは、背景の反射率の波長微分値において、背景の種類に拘わらずに同じ値とみなせるものを、基準値として求めておく手順である。
１つは、背景の反射率が基準値をとると仮定して、任意性のある解として算出された混合比及び独立成分を、それぞれ１つの解に決定する手順とである。
【００１６】
【発明の実施の形態】
以下、図面を用いて本発明の実施の形態を説明する。
【００１７】
＜本実施形態の構成＞
図１は、本実施形態の画像解析装置のブロック図である。本実施形態は、請求項１〜請求項３に対応する。画像解析装置１２は、データ取得部１６と、分析部２０と、基準値設定部２４と、任意性消去部２８と、解析結果表示部３２とで構成されている。
【００１８】
データ取得部１６は、画像データにおけるミクセル領域（画素がミクセルである画像領域）のデータ（即ち、複数の混合スペクトルで与えられるデータ）を、分析部２０に入力する（後述するステップＳ１参照）。
分析部２０は、入力された画像データにＩＣＡを施して、対象物の占有率を任意性のある解として算出する（後述するステップＳ２参照）。
【００１９】
基準値設定部２４は、背景の反射率の波長微分において安定した値を、基準値として求めておく（後述するステップＳ３参照）。
任意性消去部２８は、背景の反射率が前記した基準値をとると仮定して、対象物の占有率を、１つの解に決定する（後述するステップＳ４、Ｓ５参照）。
解析結果表示部３２は、任意性消去部２８の解析結果を、不図示の記録媒体に記録すると共に液晶モニタに表示する（後述するステップＳ６参照）。
【００２０】
＜本実施形態の解析対象モデル＞
次に、本実施形態の解析方法の説明に先立って、本実施形態の画像解析装置１２の解析対象となる周期構造の観測モデルについて説明する。
図２は、周期構造を、円形の受容野を持つセンサで観測する際のモデルの一例を示している。図において、ｘは、円形窓（ピクセル）の中心位置であり、円形窓は、互いにその外縁で接するように（重ならないように）移動する。
【００２１】
本実施形態では、観測対象領域には対象物と背景のみが存在し、同一観測対象領域内では対象物と背景のパターンは規則的であると仮定する。この場合、ピクセル位置xでの対象物の占有率をηV(x)、ピクセル位置xでの背景の占有率をηS(x)とすれば、任意のピクセル位置x1において、次式が成り立つ。
ηV(x1)＋ηS(x1)＝１  ・・・（２）
（２）式において、占有率ηV(x)とηS(x)との差が大きい場合（対象物、または背景が観測窓内の殆どを占有している場合）、ピュアスペクトルに近くなる。この場合、混合スペクトルを分離する必要がなくなる。従って、本実施形態では、どちらの占有率もある程度の大きさがある場合を想定する。
【００２２】
図３は、ピクセル位置をｘ軸とし、対象物及び背景の占有率をｙ軸としたグラフの一例である。また、図４は、図３から得られる占有率のヒストグラムの一例である。
【００２３】
図４が横軸方向に非対称であることから、周期的構造を有する対象領域上を低解像度のセンサでスキャンして得られる信号のヒストグラムは、非ガウス性であるといえる。このように構造物よりも大きな受容野を持つセンサでスキャンする代わりに、低解像度の撮像系を用いて空間的に連続したデータを取得しても、同様のヒストグラムが得られる。そして、このような非ガウス性の信号は、ＩＣＡを用いて、他のガウス性ノイズから高精度で分離できる。
【００２４】
＜本実施形態の解析方法＞
図５は、本実施形態の画像解析装置１２の解析手順を示す流れ図である。以下、図に示すステップ番号に従って、画像解析装置１２の解析手順を説明する。
【００２５】
［ステップＳ１］
ユーザは、画像解析装置１２に、画像データを入力する。データ取得部１６は、入力された画像データにおいて、ミクセルからなり、対象物及び背景の周期的な空間分布を示す画像領域があるか否かを判定する。あると判定された場合、このミクセル領域のデータを、分析部２０に入力する。
なお、データ取得部１６に機能に関しては、画像データにおけるユーザが選択した画像領域のデータを、分析部２０に入力するようにしてもよい。或いは、入力された画像データをそのまま、分析部２０に入力するようにしてもよい。
【００２６】
［ステップＳ２］
分析部２０は、以下に述べる理論に基づいて、入力されたミクセル領域のデータにＩＣＡを施し、対象物の占有率を任意性のある解として求める。
まず、ピクセル位置x、波長λでの観測信号値をI(λ,x)とすると、観測される反射光の強度は、次式で表される。
I(λ,x)＝I0(λ)×｛RV(λ)×ηV(x)＋RS(λ)×ηS(x)｝  ・・・（３）
上式において、I0(λ)は入射光スペクトル強度、RV(λ)は対象物の反射率、RS(λ)は背景の反射率である。
【００２７】
十分な解像度を得られない画像観測分野において、ＩＣＡを適用するためには、ＩＣＡの新たなモデルを構成する必要がある。そこで、本実施形態では、以下の４つに示すように、従来のＩＣＡでの各パラメータを、ミクセルを含む画像データの各パラメータに置き換えて考える。
【００２８】
第１に、従来のＩＣＡでの複数の観測点の位置を、波長（λ）とする。
第２に、従来のＩＣＡでのサンプリング点を、画像データ内での画素の位置、即ち、観測窓の位置（ｘ）とする。
第３に、従来のＩＣＡでの音の原信号を、位置（ｘ）の観測窓内における、対象物の占有率ηV(x)及び背景の占有率ηS(x)とする。
第４に、従来のＩＣＡでの混合比を、各波長における、対象物の反射率RV(λ)と背景の反射率RS(λ)との差とする。
【００２９】
なお、ＩＣＡを適用するためには、以下の３つの条件を満たす必要がある。
第１に、混合信号は、原信号（占有率ηV(x)、ηS(x)）の線形和である。
第２に、原信号の確率密度分布は、非ガウス性である。
第３に、各構成要素である原信号は、それぞれ独立である。
図４で説明したように、第２の条件は満たされている。また、本実施形態では、（２）式が満たされるモデルを想定しているため、第１及び第３の条件も満たされている。
【００３０】
ここで、（２）式を（３）式に代入すると、次式が得られる。
I(λ,x)＝I0(λ)×［｛RV(λ)−RS(λ)｝×ηV(x)＋RS(λ)］・・・（４）
上式から明らかなように、本実施形態の解析対象モデルでは、１つの構成要素（この例では対象物）のみを推定することと等価になる。
また、実際の観測条件下では、撮像素子の熱雑音、撮像感度のばらつき、大気ノイズ等のノイズ成分を考慮する必要がある。従って、次式のように、（４）式にノイズ項を加算する必要がある。なお、ノイズ成分は通常、ガウス性である。
【００３１】
2004184391.tif 000005
上式において、Rn(λ)はノイズの振幅であり、ηn(x)はノイズの周期パターンである。（５）式における、これらRn(λ)とηn(x)との積は、請求項記載のノイズ成分、N(λ，x)に対応する。また、実際には占有率（被覆率）にも各対象物サイズの個体差による揺らぎが存在するため、次式のように揺らぎ成分を考慮する必要がある。
【００３２】
ηV(x)＝ηVtrue(x)＋ΔηV(x)  ・・・（６）
また、揺らぎ成分を除いた対象物の占有率ηVtrue(x)を、交流成分ηVAC(x)と直流成分ηVDCに分け、次式のように表す。なお、直流成分ηVDCは、ピクセル位置ｘに依存しない。
ηV(x)＝ηVAC(x)＋ηVDC＋ΔηV(x)  ・・・（７）
従って、ノイズ成分や被覆率の揺らぎを考慮した場合、（７）式を（５）式に代入することで、観測される混合スペクトルは次式のように表される。
【００３３】
2004184391.tif 000006
次に、前記したＩＣＡの基本となる（１）式に基づいて、（８）式を次式のように定義する。なお、本発明では、ＩＣＡで算出される任意性を有する解を一意に定めるために、z0をピクセル位置xに対する定数項成分（直流成分）として扱う。
【００３４】
2004184391.tif 000007
ＩＣＡを適用する場合、周期性を示す信号（xに依存する交流成分）は独立成分（原信号）sとして見なされ、周期性を持たない信号（xに依存しない直流成分）は定数項成分z0と見なされる。即ち、（９）式は、以下の２式のように、周期性を示す独立成分sと定数項成分z0とに分離可能と考えられる。
【００３５】
2004184391.tif 000008
ここで、光源の状態により観測値が異なると評価が困難である。そこで、観測信号値を、入射光強度に対する比率として求めて、求めた比率をI(λ,x)とする。このようにすれば、I0(λ)＝１とみなせる（但し、以下の式でも便宜上、I0(λ)を残して記載する）。次に、（１０）式に対して、次式のようにＩＣＡを適用する。
【００３６】
2004184391.tif 000009
上式において、ＩＣＡにより求まるAVAC、AnAC、sVAC、snAC は、それぞれ以下の項に対応する。
【００３７】
混合比の行列Aにおいて交流成分に対応するAVAC は、
I0(λ)×｛RV(λ)−RS(λ)｝ に対応する。
独立成分の行列sにおいて交流成分に対応するsVAC は、
｛ηVAC(x)＋ΔηV(x)｝に対応する。
混合比の行列Aにおいてノイズ成分に対応するAnAC は、Rn(λ)に対応する。
独立成分の行列sにおいてノイズ成分に対応するsnAC は、ηn(x)に対応する。
【００３８】
ＩＣＡでは、推定される複数の原信号の絶対強度と、複数の原信号の混合比を表す混合行列の各要素の大きさのとり方に任意性が残されている。即ち、（１２）式を満たすという条件において、交流成分に対応する混合比と独立成分とは、ＩＣＡにより算出された１つの解（AVAC、sVAC）以外に、幾つかの解を持つことができる。この任意性は、ＩＣＡでは白色化処理のために値を正規化（平均ゼロ、大きさ１に変換）することにより生じる。
【００３９】
従って、反射率RV(λ)、RS(λ)、占有率ηV(x)、ηS(x)を特定するためには、混合比AVAC及び独立成分sVACが持つ任意性を、一意に決定する必要がある。そこで、求めるべき真の混合比をnewAVAC、真の独立成分newsVACを、次式ように定義する。
newAVAC＝K×AVAC ＝I0(λ)×｛RV(λ)−RS(λ)｝・・・（１３）
newsVAC＝（１／K）×sVAC  ・・・（１４）
上式において、Kは較正係数である。そして、次のステップでは、較正係数Kを求めて、任意性を有する混合比AVAC及び独立成分sVACを、１つの解に定める方法を説明する。
【００４０】
［ステップＳ３］
較正係数Kを決定するためには、混合比AVACが一定で安定した波長λで決定する方法と、波長微分を用いて混合比AVACがより安定した波長λで決定する方法とが考えられる。本実施形態では、後者の手法を用いる。即ち、対象物の種類に依存しない波長λの微分値を特定することで任意性を消去する。
【００４１】
対象物の反射率RV(λ)は様々な変化をするが、その微分値が０となる波長はほぼ決まっている。一方、背景の反射率RS(λ)は、スペクトルが滑らかに変化する特徴がある。このため、背景の反射率RS(λ)の波長微分値は、通常一定値となる。ここで、RS(λ)の波長微分値で最も安定したものを用いてもよいが、例えば背景が土の場合、土壌種類や水分含有量等によってばらつきがある。従って、特定の波長でのRS(λ)の波長微分値のみを用いると、微分値が変動するおそれがある。
【００４２】
そこで、本実施形態では、まず、背景の反射率RS(λ)の波長微分値が一定値であり、且つ、対象物の反射率RV(λ)の波長微分値がゼロである条件が成立する波長を複数選択する。ここで複数の波長の値を利用するのは、様々な種類及び状態の背景に対して安定した値を得るためである。
次に、ミラー行列による判別分析を実施する。これは、選択した複数の波長での微分値に生じる分散を小さくし、且つ、背景の反射率RS(λ)の波長微分値の絶対値をできる限り大きくして誤差影響を軽減するためである。ここでの『誤差影響』とは、『ゼロと仮定した対象物の反射率RV(λ)の波長微分値が正確にゼロではないことによる誤差』である。また、ここでの『絶対値をできる限り大きく』する方法としては、元の行列に対して、座標の原点を中心として対称位置にくるミラー行列を用いることが有効と考えられる。
【００４３】
なお、背景の反射率RS(λ)の波長微分値を求める手段としては、発光波長の僅かに異なる複数のレーザ光を照射し、それぞれの波長で得られる画像間の差分をとる方法が考えられる。或いは、波長選択性のない光源を用い、撮影レンズに装着する光学フィルタを切り替えて複数の画像を撮影し、その差分から波長微分値を算出してもよい。このとき、複数の光学フィルタは、透過波長域特性が僅かに異なるものを用いる。
【００４４】
以上の理論に基づいて、基準値設定部２４は、背景の反射率の波長微分において安定した値（後述するｎｅｗＲｅｆSJ’、及び請求項記載の基準値に対応）を、以下のように求めておく。
まず、対象物の反射率RV(λ)の微分値が０となる波長λi(i＝１,２,・・・,p) における、特定の１つの背景の反射率をＲｅｆSJと定義する。ここで、Ｊは背景の種類であり、Ｊ＝１、２、・・・ｍとする。ＲｅｆSJは、次式で表されるように、p行１列である。
【００４５】
ＲｅｆSJ＝｛RefSJ(λ1), RefSJ(λ2),・・・, RefSJ(λp)｝T・・・（１５）
また、特定の１つの背景の反射率ＲｅｆSJの波長微分値をＲｅｆSJ’と定義し、これを次式で表す。
ＲｅｆSJ’＝｛RefSJ'(λ1), RRefSJ'(λ2),・・・, RefSJ'(λp)｝T・・・（１６）
次に、上式に基づいて、背景の反射率の（p×ｍ）個の微分値分布をClass1と定義する。Class1は、次式で表されるように、p行ｍ列である。
【００４６】
Class1＝｛ＲｅｆS1’, ＲｅｆS2’,・・・, ＲｅｆSm’｝・・・（１７）
また、Class1のミラー行列を、Class2と定義する。Class2は、次式で表される。
2004184391.tif 000010
これらClass1、Class2に対して、統計学で用いられている判別分析を適用する。なお、判別分析は、２つのグループＸ、Ｙに予め分けられたデータがあるときに、新たに入ってきたデータがどちらのグループに属するかを判定する方法である。判別分析では、グループＸ、Ｙは、それぞれ、軸変換の行列を乗じられて、別の座標系で表されたグループＸ’、Ｙ’に写像される。判別分析の効果として、グループＸ’、Ｙ’は、グループＸ、Ｙよりもグループ内の分散が小さくなる。また、グループＸ’、Ｙ’間の距離は、グループＸ、Ｙ間の距離よりも大きくなる。
【００４７】
ここで、Class1、Class2に対して軸変換の行列を乗じて写像した後の座標系の軸を、Λとする。なお、写像前のClass1、Class2の座標系の軸は、波長λである。また、写像後のClass1をClass1’とし、写像後のClass2をClass2’とする。行列Class1’の行数及び列数は、写像前と同じである。同様に、行列Class2’の行数及び列数は、写像前と同じである。そして、Class1’の各列ベクトルを、次式のように定義する。
【００４８】
2004184391.tif 000011
上式右辺における１（または２、３、・・・）列目の全要素は、次式において、Ｊ＝１（または２、３、・・・）  としたものと定義する。
ｎｅｗＲｅｆSJ’＝｛newRefSJ'(Λ1),・・・, newRefSJ'(Λp)｝T・・・（２０）
判別分析により得られるClass1’の各要素ｎｅｗＲｅｆSJ’ は、背景の反射率の微分値において、背景の種類に影響を受けにくく安定したものとなる（背景の種類に拘わらずほぼ同じ値となる）。
【００４９】
［ステップＳ４］  任意性消去部２８は、例えば、Ｊ＝１のときのｎｅｗＲｅｆSJであるｎｅｗＲｅｆS1’を用いて、以下のように較正係数Kを求める。なお、Ｊ＝１に限らず、他の値（２〜ｍ）でもよい。
まず、（１３）式において、RS(λ)がｎｅｗＲｅｆS1’に相当し、R(λ)が後述するｎｅｗＲｅｆV1’に相当すると見なせば、次式が成り立つ。
【００５０】
2004184391.tif 000012
前記したように、対象物の反射率RV(λ)の微分値がゼロであることを仮定しているので、ｎｅｗＲｅｆV1’はゼロとみなしてよい。なお、ｎｅｗＲｅｆV1’を具体的に求めるならば、まず、ＲｅｆSJ’と同様に、対象物の反射率の微分値ベクトルを次式のように定義する（次式で、Ｊ＝１、２、・・・ｍ）。
【００５１】
ＲｅｆVJ’＝｛RefVJ(λ1), ・・・, RefVJ(λp)｝T・・・（２２）
次に、ＲｅｆVJ’を１列分の要素とするp行ｍ列の行列をGroup1とし、そのミラー行列をGroup2とする。そして、これら Group1、Group2 に対して判別分析を施して同様に写像する。ｎｅｗＲｅｆV1’ は、写像後のGroup1における１列目の要素に対応する。
（２１）式において、AVACと、I0(λ)と、ｎｅｗＲｅｆV1’と、ｎｅｗＲｅｆS1’とは既に求めたから、較正係数Kを算出できる。そして、任意性消去部２８は、求めた較正係数Kを（１３）、（１４）式に代入して、混合比newAVAC及び独立成分newsVACを算出する。
【００５２】
［ステップＳ５］
任意性消去部２８は、（９）式において、スキャン方向ｘについてI(λ,x)の平均値を、mean｛I(λ,x)｝として求める。
これにより、次式のように、交流成分ηVAC(x)は消去され、揺らぎ成分ΔηV(x)及びガウス性ノイズ成分Rn(λ)×ηn(x)はゼロに収束する。
【００５３】
mean｛I(λ,x)｝＝I0(λ)×［｛RV(λ)−RS(λ)｝×ηVDC＋RS(λ)］・・・（２３）
そこで、（２３）式に、（１３）式を代入する。
2004184391.tif 000013
任意性消去部２８は、上式において、較正係数Kを決定したときと同様に、RS(λ)の１つの値として例えばｎｅｗＲｅｆS1’を代入して、直流成分ηVDCを求める。そして、任意性消去部２８は、以下の５つを（９）式に代入して、背景の反射率RS(λ)を求める。
【００５４】
第１に、本ステップで求めた直流成分ηVDC。
第２に、ステップＳ４で求めた混合比newAVAC を、
I0(λ)×｛RV(λ)−RS(λ)｝ として代入する。
第３に、ステップＳ４で求めたsVACを、｛ηVAC(x)＋ΔηV(x)｝として代入する。
第４に、ステップＳ２において（１２）式にＩＣＡを適用したときに求めた
AnACを、Rn(λ)に代入する。
第５に、ステップＳ２において（１２）式にＩＣＡを適用したときに求めた
snACを、ηn(x)に代入する。
これにより、背景の反射率RS(λ)が求まる。従って、newAVACとして求めた
I0(λ)×｛RV(λ)−RS(λ)｝を用いて、対象物の反射率RV(λ)も求まる。
【００５５】
［ステップＳ６］  解析結果表示部３２は、以上の解析により得られた較正係数K、混合比newAVAC、独立成分newsVAC、直流成分ηVDC、背景の反射率RS(λ)、対象物の反射率RV(λ)を、記録媒体に記録すると共に液晶モニタに表示する。
【００５６】
図６及び図７は、信号強度よりも強いガウス性ノイズが混入された、図２のようにある程度のばらつきを有する周期構造の観測画像に、上述の解析処理を施したシミュレーション結果である。図６はノイズから分離された混合スペクトルを示しており、図７は分離されたガウス性ノイズ成分を示している。なお、対象物として植物、背景として土壌を想定している。
【００５７】
＜本実施形態の効果＞
本実施形態では、対象物の反射率RV(λ)と背景の反射率RS(λ)との差を混合比AVACとすると共に、対象物の占有率η(x)を独立成分sVACとする。次に、ＩＣＡにより、混合信号I(λ,x)から、混合比AVAC及び独立成分sVACを任意性のある解として算出する（ステップＳ２）。
【００５８】
次に、背景の反射率RS(λ)の微分値が一定値であり、且つ、対象物の反射率RV(λ)の微分値がゼロである条件が成立する波長を複数選択する。そして、選択した波長での背景の反射率RS(λ)の微分値からなる行列と、この行列のミラー行列に判別分析を実施する。これにより、背景の種類に影響を受けにくく安定した、背景の反射率の微分値ｎｅｗＲｅｆSJ’を求める（ステップＳ３）。
【００５９】
そして、ｎｅｗＲｅｆSJを用いて、混合比及び独立成分の任意性を一意に定める較正係数Kを求める（ステップＳ４）。従って、画像データにおける、ミクセルからなり、対象物の周期的な空間分布を示す画像領域のデータから、対象物の反射率RV(λ)、及び対象物の占有率ηV(x)を算出できる。
【００６０】
また、ＩＣＡを用いるので、画像データに含まれるガウス性ノイズ成分が大きくても、ノイズ成分を分離して、周期的な空間分布を示す対象物の占有率を高精度で推定できる。なぜなら、空間に周期的に位置する対象物の占有率分布は、非ガウス性であるからである。
本実施形態の手法は、例えば以下に挙げるような、周期的分布を持つ対象領域に適用できる。内視鏡を用いた体内組織像の観測、光学顕微鏡の分解能限界を超えたナノ周期構造の観測、周期信号に変調を与える形で書き込まれた種々の微細記憶媒体からの信号の読み出し、遺伝子情報の解析、都市部の環境モニタリングなどである。
【００６１】
＜本実施形態の補足事項＞
なお、上述したステップＳ１〜Ｓ６の処理をプログラムコード化して画像解析プログラムを作成し、コンピュータ上で上述した画像解析を実施してもよい。この場合、請求項４に対応し、上述した実施形態と同様の効果を得ることができる。また、ステップＳ１〜Ｓ６で説明した画像解析方法は、請求項５にも対応する。
【００６２】
【発明の効果】
本発明では、周期的な空間分布を有する対象物の反射率及び背景の反射率の混合信号から、ＩＣＡにより、対象物の占有率を任意性のある解として求める。そして、背景の反射率が、背景の反射率の波長微分において安定した値をとると仮定する。これにより、前述の任意性を１つの解に定める。従って、前述の混合信号から、対象物の占有率を算出できる。
【図面の簡単な説明】
【図１】本発明の画像解析装置のブロック図である。
【図２】周期構造を、円形の受容野を持つセンサで観測する際のモデルの一例を示す説明図である。
【図３】ピクセル位置をｘ軸とし、対象物及び背景の占有率をｙ軸としたグラフの一例である。
【図４】図３から得られる占有率のヒストグラムの一例である。
【図５】本発明の画像解析装置の解析方法を示す流れ図である。
【図６】ガウス性ノイズを含む周期構造の観測画像に本発明の解析処理を施して得られた、ノイズから分離された混合スペクトルの波形図である。
【図７】図６のシミュレーション結果における、分離されたガウス性ノイズ成分の波形図である。
【符号の説明】
１２  画像解析装置
１６  データ取得部
２０  分析部
２４  基準値設定部
２８  任意性消去部
３２  解析結果表示部
【図１】
2004184391.tif 000015
【図２】
2004184391.tif 000016
【図３】
2004184391.tif 000017
【図４】
2004184391.tif 000018
【図５】
2004184391.tif 000019
【図６】
2004184391.tif 000020
【図７】
2004184391.tif 000021