森林の樹冠の評価方法及びその樹冠評価プログラム
【課題】 森林の高分解能の衛星写真画像から樹冠形状を検出することができ，その樹冠の樹種を正確に判定することができる森林の樹冠評価方法及びその樹冠評価プログラムを提供する。
【解決手段】 森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠評価方法において，前記画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行い，その平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，樹冠の形状を求める。更に，樹冠の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する。
【解決手段】 森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠評価方法において，前記画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行い，その平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，樹冠の形状を求める。更に，樹冠の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は，森林の樹冠の評価方法及びその樹冠評価プログラムに関し，特に，衛星写真画像など森林を上空から撮影した高分解能の画像データに対して所定の画像処理を行って森林の樹冠の位置，形状（大きさ）を検出し樹種を判定する樹冠評価方法及びその樹冠評価プログラムに関する。
【背景技術】
【０００２】
森林地帯の衛星写真画像データに画像処理を加えることで森林の樹冠を評価することが注目されている。樹冠とは，森林を構成する樹木を上部から観察した時に，樹木の枝葉で構成された丸みを持った形状をいい，樹種とは樹木の種類をいう。森林を構成する樹冠の位置，形状，大きさの検出に加えてその樹木の種類（樹種）を判定することで，例えば二酸化炭素の排出権の評価や，森林の評価を行うことができる。
【０００３】
従来から，人工衛星や航空機から撮影した森林のカラー写真画像を画像処理して，森林の樹冠形状画像を求め，その樹冠形状画像と，実際に調査した樹冠に対応するカラー画像とを比較して，樹冠の樹種を判定する森林評価方法が提案されている。例えば，以下の特許文献１に記載されている。
【０００４】
この森林評価方法によれば，人工衛星や航空機で撮影した森林のカラーデジタル写真画像を白黒デジタル画像に変換し，平滑化処理し，樹冠の境界線の強調処理をし，ウォーターシェッドアルゴリズムにより樹冠形状を求め，樹冠形状画像内の画素のスペクトル分析により樹種を判定する。
【０００５】
また，人工衛星から撮影した森林地帯の写真画像から，樹木の活性度の指標である正規化植生指数（ＮＤＶＩ：Normalized Difference Vegetation Index）を算出する森林の評価方法や，森林の写真画像の濃淡パターンに対して，その一様性，方向性，コントラストなどを濃度共起行列を用いて求め，その特徴量から森林の樹種を判定する方法や，ヘリコプタから撮影した森林の写真画像に対して樹冠形状を近似円と仮定してその近似円を森林の写真画像に適合させて樹冠を検出する方法なども提案されている。例えば，以下の非特許文献１に記載されるとおりである。
【０００６】
更に，高分解能衛星，例えば地球観測衛星ＩＫＯＮＯＳ，の画像を利用した植生解析が日本スペースイメージング社から提案されている。これによれば，高分解能衛星画像を前処理し，領域分割して樹冠形状を求め，各分割領域のスペクトル解析を行って，その波長スペクトルを基にして樹種を分類することが提案されている。但し，領域分割の方法については公表されておらず，スペクトル解析の具体的方法も公表されていない。
【０００７】
そして，衛星ＩＫＯＮＯＳからの高分解能の森林の画像を現地調査で同定させた樹冠形状に重ねて，その樹冠の画素の平均値と標準偏差を計算し，その樹冠のスペクトル分析により樹種を分類することが提案されている。例えば，以下の非特許文献２に記載されるとおりである。
【特許文献１】特開２００１−３５７３８０号公報
【非特許文献１】「衛星リモートセンシングによる植生モニタリング−有峰地域の森林活性度調査−平成１２−１４年度共同研究成果報告書」平成１５年３月，金沢大学工学部村本研究室，北陸電力株式会社立地環境部
【非特許文献２】「高分解能ＩＫＯＮＯＳ衛星による針広混交林の樹種分類」加藤正人，森林航測 大１９８号 ６〜９ ２００２．１１
【発明の開示】
【発明が解決しようとする課題】
【０００８】
非特許文献１に記載されている衛星画像（例えばＬａｎｄｓａｔ）のような画素当たり３０ｍと低分解能の場合は，画素内に複数の樹冠が含まれ，そもそも樹冠形状を検出することができない。したがって，森林全体の評価をするに止まる。また，同文献に記載されているような航空機写真画像の場合は，樹冠検出は近似円による方法が提案されているにすぎず，また写真の濃淡パターンによる樹種判定では，樹木が密集した森林全体の樹種判定をするに止まっている。
【０００９】
更に，特許文献１に記載された衛星画像に平滑化処理等を行ってウォーターシェッドアルゴリズムにより樹冠形状を求める方法によると，輝度変化が最も小さい部分を中心にして，輝度の変化に沿って四方に領域を拡大させ，隣接する拡大領域が衝突する箇所を樹冠の境界線と定めるものであるが，樹冠の境界部分の輝度値変化を残す程度の平滑化処理だけでは、高分解能の衛星画像データの樹冠の輝度は単純な峰の形状にならないので，樹冠を正確に検出できない。また，樹種の判定に樹冠の画像のスペクトル解析を利用する方法が記載されているが，樹冠画像が複数の画素からなる高分解能の画像データの場合に波長スペクトルの相関性を利用した方法を適用すると，各画素が同じ波長スペクトルを有するとは限らず，正確に樹種を判定することは困難である。
【００１０】
また，非特許文献２では，高分解能の衛星写真画像を利用しているものの，樹冠形状の検出は，実地調査画像との重ね合わせによるものであり，実地調査自体が高コストであり現実的な方法とはいえない。また，スペクトル解析による樹種判定は，前述のとおり正確ではない。
【００１１】
そして，日本スペースイメージング社の提案では，前述のとおり樹冠形状の検出方法やスペクトル解析による樹種判定の具体的方法については記載されていない。
【００１２】
そこで，本発明の目的は，森林の高分解能の衛星写真画像から樹冠形状を正確に検出することができ，その樹冠の樹種を正確に判定することができる森林の樹冠評価方法及びその樹冠評価プログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記の目的を達成するために，本発明の第１の側面によれば，森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠評価方法において，前記画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行い，その平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，樹冠の形状を求める。更に，樹冠の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する。
【００１４】
本発明の第１の側面において，好ましい実施例では，平坦化処理が，前記画像データの輝度値を所定値だけ減算し元の輝度値を超えない範囲で近傍画素の最大輝度値に置き換える膨張処理と，画像データの輝度値を所定値だけ加算し元の輝度値を下回らない範囲で近傍画素の最小輝度値に置き換える縮退処理とを有する。この実施例によれば，コンピュータソフトウエアを利用した画像処理により膨張処理と縮退処理を簡単に行うことができ，それにより峰と谷が平坦化された画像データについて領域分割処理を簡単に行うことができる。
【００１５】
更に，好ましい実施例では，領域分割処理は，前記平坦化処理された画像データの微分値を求める微分値生成処理と，微分値ゼロや極小値の位置を中心に隣接する画素領域を拡張する領域拡張処理と，極小値の位置から拡張された領域を微分値ゼロの位置から拡張された領域に併合する併合処理とを有する。微分値を利用することにより，コンピュータソフトウエアを利用した画像処理で領域分割処理を簡単に行うことができる。
【００１６】
本発明の第１の側面において，好ましい実施例によれば，上空から撮影した画像データは，単一画素当たり数１０センチメートル〜数メートル程度の分解能であり，単一樹冠内に複数行複数列の画素が含まれ且つ単一画素内に複数の枝葉が含まれる程度の分解能である。かかる分解能を有する画像データの場合，樹冠画像には樹種に応じた固有の模様（テクスチャ）が生成されるので，スペクトル分析でなくテクスチャ特徴量による樹種判定を容易に行うことができる。
【００１７】
更に，好ましい実施例では，前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点の輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・・，ｎ−１）を要素とする同時生起行列から求められる第１のテクスチャ特徴量，画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含む第２のテクスチャ特徴量，画像内で一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められる第３のテクスチャ特徴量，画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められる第４のテクスチャ特徴量のいずれかを含む。これらのテクスチャ特徴量を利用することにより，適切に樹種を判定することができる。
【００１８】
上記の目的を達成するために，本発明の第２の側面は，第１の側面の評価方法をコンピュータに実行させる樹冠評価プログラムである。
【発明の効果】
【００１９】
上記本発明の第１の側面によれば，樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行うことにより，その後の領域分割処理において，ウォーターシェッドアルゴリズム（流域検出方法）による簡単な画像処理により樹冠形状の検出を正確に行うことができる。また，樹冠の画像のテクスチャ特徴量により樹種を判定するので，樹冠内の複数の画素のスペクトル分析よりもより正確に樹種を判定することができる。
【発明を実施するための最良の形態】
【００２０】
以下，図面にしたがって本発明の実施の形態について説明する。但し，本発明の技術的範囲はこれらの実施の形態に限定されず，特許請求の範囲に記載された事項とその均等物まで及ぶものである。
【００２１】
図１は，本実施の形態における人工衛星など上空から撮影した画像と樹冠との関係を示す図である。この図により本実施の形態が対象とする画像データの分解能を説明する。森林を上空から撮影すると，樹木の枝葉で構成された丸みを持った樹冠１０，１２，１４，１６が観察される。衛星ＩＫＯＮＯＳなどの高分解能衛星の画像データの場合は，その分解能が１画素当たり１メートル，または６０センチメートルであり，樹冠の直径が数メートルから数十メートルとすると，図示されるとおり各樹冠は複数行複数列の画素を含むことになる。図１の例では，実線の四角形が樹冠の画素に対応し，破線の四角形が樹冠以外の画素に対応する。また，１メートル／画素程度の分解能では，１画素内に樹木の複数の枝葉が含まれる。したがって，樹木の枝葉を複数の画素で表現するほどの高い分解能ではない。
【００２２】
このように，衛星ＩＫＯＮＯＳによる画像データは，樹冠の形状を目視判読できる程度に高い分解能を持つので，かかる画像データをコンピュータを利用して画像処理することによりその樹冠形状を検出することができる。そして，検出された樹冠形状の樹種を自動判定できれば，森林の評価方法としてその価値を高めることができる。近年，地球温暖化に伴う二酸化炭素の排出権ビジネスの始まりに際して，植林により創出される排出権を効率的且つ正確に評価するツールが求められているところであり，かかる衛星画像データを利用した自動森林評価方法は，排出権取引に欠かせないものである。特に，衛星画像データから樹木１本１本の形状を検出し，その樹種を判定できれば，正規化植生指数（ＮＤＶＩ：Normalized Difference Vegetation Index）による植生活性度（樹木の葉緑素の度合い）などの算出と組み合わせるなどして，より正確に二酸化炭素の吸収量を評価することができる。
【００２３】
図２は，本実施の形態における森林の樹冠の評価方法についてのフローチャート図である。この評価方法は，森林を上空から撮影した画像データをコンピュータプログラムにより画像処理を行うことで自動的に行うことができる。したがって，この評価方法は，汎用コンピュータに当該画像処理プログラムをインストールした評価コンピュータにより行うことができる。
【００２４】
また，図３は，上記評価方法における森林の樹冠形状の検出方法のフローチャート図である。これらのフローチャートにしたがい，且つ具体的画像例を参照しながら，森林の樹冠の評価方法を説明する。
【００２５】
図２には，森林樹冠の評価方法の概略が示されている。まず，衛星画像データを入力する（Ｓ１０）。この衛星画像データは，青，緑，赤，近赤外（Ｂ，Ｇ，Ｒ，ＮＩＲ）の各波長域の地表反射データからなるマルチスペクトル画像データである。図４，５に，これらの画像の一例を示す。図４には，青と緑の衛星画像が示される。また，図５には，赤と近赤外の衛星画像が示される。各画像データは，画素に対応する輝度値であり，例えば８ビット，２５６階調の輝度値データである。これらの衛星画像から目視判読されるとおり，複数の画素を含む樹冠が検出可能な程度の分解能を有している。
【００２６】
次に，４つの画像データについて主成分分析して，第１主成分の画像を生成する（Ｓ１２）。この工程については，後に詳述するが，要すれば，図４，５に示した４つの画像データを４次元の輝度値と考えて，それら４次元の輝度値から最も特徴の現れる軸（特徴の分散が最大になる軸）を求め，４つの画像データを線形変換して当該軸上の輝度値（第１主成分）を求めるものである。
【００２７】
そして，第１主成分の画像データに対して膨張処理と縮退処理を行い，画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更することなく，峰と谷の部分を削除して平坦にする（Ｓ１４）。この平坦化処理を行う前の第1主成分の画像データの輝度値の空間変化では，樹冠の境界部分の輝度値が大きく増減すると共に，樹冠内においても輝度値が増減する峰と谷が複数存在する。このような峰や谷は，ウォーターシェッドアルゴリズムによる領域分割処理を行うとそれぞれ個別の領域となり、一つの樹冠が細分化される。この細分化を防ぎ正確な樹冠形状を検出するためには，平坦化処理により峰や谷を平坦化する処理が有効である。
【００２８】
平坦化処理された画像データの輝度値の空間変化に対して，ウォーターシェッドアルゴリズム（流域検出方法）を適用して，樹冠の領域分割処理を行い，樹冠の位置や形状（領域）を求める（Ｓ１６）。ウォーターシェッドアルゴリズム自体は公知の手法であるが，そのうち，本実施の形態では，画像データの輝度値の空間変化について微分値を求め，微分値ゼロや微分値の極小値の位置を開始点として，周囲の隣接する画素領域を拡張する領域拡張処理を行う。平坦化された画像データの微分値を利用することで，ウォーターシェッドアルゴリズムによる領域分割処理をコンピュータプログラムにより容易に行うことができる。
【００２９】
樹冠形状が検出された後，その樹冠の画像データ，ここでは第１主成分の画像データ，の空間的模様であるテクスチャ特徴量を抽出し，既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較し，樹種を判定する（Ｓ１８）。本実施の形態では，衛星ＩＫＯＮＯＳの画像の分解能では樹冠の画像データが複数の画素により構成されることを考慮して，従来のスペクトル分析ではなく，テクスチャ特徴量を利用して樹種判定を行う。このテクスチャ特徴量は，例えば，同時生起行列から求められる局所一様性などの特徴量が有効である。
【００３０】
そして，最後に評価対象の森林における樹冠の位置，形状（大きさ），樹種を出力する（Ｓ２０）。この出力データと前述のＮＤＶＩ値を利用して，森林の排出権の評価を行うことができる。
【００３１】
［樹冠の検出］
次に，図３のフローチャートにしたがって，森林の樹冠の検出方法について詳述する。樹冠検出は，図２の工程Ｓ１０〜Ｓ１６により行われ，図３にはそれらの構成を詳述しており，同じ工程には同じ引用番号を与えている。まず，衛星画像データの入力（Ｓ１０）は，図２と同じである。次に，第１主成分の画像生成工程Ｓ１２が行われる。
【００３２】
図６は，第１主成分の画像生成を説明するための図である。第１主成分の画像を生成する主成分分析については，例えば，「新編，画像解析ハンドブック，高木幹雄・下田陽久監修」に詳しく説明されており，それによると，主成分分析とは，多変量の測定値（４色の輝度値）が与えられたとき，線形変換により変量間の相関をなくし，より少ない変量により測定対象の特徴を記述しようとする変換である。図６により簡単に説明すると，仮に２種の変量1ｘ，2ｘが測定された時，二次元座標（1ｘ軸と2ｘ軸）上には，黒点の変量２２がプロットされる。そして，これら変量の分散が最大になる破線の軸が主成分軸２０になる。
【００３３】
本実施の形態の４色の画像データの輝度値に適用して説明すると，図６の式（１）は，式（３）の条件の下で変量Ｚの分散が最大になるように，式（２）の係数ａを定めた場合の，第１主成分である変量Ｚを示す。この変量Ｚは，４色の輝度値である変量1ｘ，2ｘ，3ｘ，4ｘを係数ａ1,a2,a3,a4との一次結合で表したものである。
【００３４】
この第１主成分Ｚは，４色の輝度値について共分散行列の固有値方程式を解き，その最大固有値に対応する固有ベクトルを第１主成分ベクトルａとし，式（１）のとおり４つの輝度値1ｘ，2ｘ，3ｘ，4ｘに第１主成分ベクトルａを乗算して求める。図７は，図４，５の４つの画像データから求めた第１主成分の画像である。各画素の輝度値は，その分散が最大になっており，樹冠形状や樹冠の模様（テクスチャ特徴量）を求めるのに適した値に変換されている。
【００３５】
図３に戻り，次に，第１主成分の画像に対して，平滑化処理を行う（Ｓ１３）。この平滑化処理は，従来から知られているガウシアンフィルタを利用した処理である。ガウシアンフィルタの一例が図８に示されている。この例では，ガウシアンフィルタＧＦは３×３のフィルタであり，平滑化処理では，画像データの各画素の輝度値とその周囲８近傍の画素の輝度値に対して，このガウシアンフィルタＧＦの係数をそれぞれ乗算し，その合計値が中心の画素の輝度値に置き換えられる。これにより，第１主成分の画像データの空間変化における高周波成分が除去され平滑化される。このガウシアンフィルタＧＦは，８近傍の画素の輝度値を対応する係数に応じて加算するものであり，そのサイズを３×３と小さくすることで，樹冠を検出するための境界部分の輝度値変化を残しつつ，平滑化することができる。
【００３６】
図８は，図６の第１主成分の画像を上記ガウシアンフィルタにより平滑化処理した画像を示す図である。画像データの空間変化が平滑化されて，高周波成分が除去されているが，樹冠を検出するための境界部分の輝度値の大きな変化（低周波成分）は残されたままである。
【００３７】
次に，平坦化処理Ｓ１４のために，膨張処理により輝度値の空間変化の峰の部分を除去して平坦化し（Ｓ１４Ａ），縮退処理により輝度値の空間変化の谷の部分を除去して平坦化する（Ｓ１４Ｂ）。膨張処理と縮退処理はいずれを先に行っても良い。
【００３８】
図９は，膨張処理のフローチャート図であり，図１０は，膨張処理を説明する輝度値の空間変化のグラフ図である。また，図１１〜図１３は，図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。膨張処理は，平滑化処理された画像データの輝度値の空間変化について，その峰の部分を平坦化する処理である。図１０を参照しながら図９のフローチャートを説明すると，まず，輝度値の空間変化の峰の高さを測定し，比較的多数の峰の高さを含む高さを所定値に設定する（Ｓ２０）。図１０の実線は，平滑化処理された画像データの輝度値の空間変化を示す。峰とは極大値をとる部分であり，この例では４つの峰が存在している。この場合，峰Ｍ１とＭ２が一つの樹冠，峰Ｍ３が樹冠以外，峰Ｍ４が別の樹冠と予測される。そこで，膨張処理により，峰Ｍ１とＭ２を平坦化して一つにし，峰Ｍ３はなくし，峰４は上部の一部分のみを平坦化する。そのために，極大値の位置にある峰の高さを近傍の谷からの高さとして検出し，それらのうち比較的多数の峰の高さを含む高さを上記所定値と設定する。ここでは，峰Ｍ１，Ｍ２，Ｍ３のうちもっと高い高さに設定される。
【００３９】
次に，各画素の輝度値を前記所定値だけ減算した画像データを生成する（Ｓ２２）。その結果，図１０中の破線のグラフ（輝度値の空間変化）が得られる。この減算した画像データの輝度値を，元の輝度値を超えない範囲で８近傍の画素の最大輝度値に置き換える処理を，全ての画素に対して行う（Ｓ２４，Ｓ２６，Ｓ２８，Ｓ３０）。つまり，注目画素の輝度値が近傍８画素の最大輝度値よりも小さい場合は（Ｓ２４のＹＥＳ），更に，最大輝度値が元の輝度値より大きくないときは（Ｓ２６のＮＯ）注目画素の輝度値を最大輝度値に置き換え（Ｓ２８），最大輝度値が元の輝度値より大きいときは（Ｓ２６のＹＥＳ）注目画素の輝度値を元の輝度値に置き換える（Ｓ３０）。上記の処理を全ての画素について行う（Ｓ３２）。この置き換え処理Ｓ２４〜Ｓ３２は所定回数繰り返すことが望ましい。
【００４０】
その結果，図１０中の一点鎖線のグラフの輝度値が得られる。つまり，実線矢印のように減算処理された破線のグラフが，破線矢印のように，近くの峰の輝度値まで（但し元の輝度値以下）上昇され，最終的に小さな峰Ｍ１とＭ２が平坦化されて一つになり，Ｍ３が平坦化され，大きな峰Ｍ４の上部一部分が平坦化される。つまり，樹冠の境界部分の急激な輝度値変化領域３０を実質的に変更せず，小さな峰や大きな峰の上部一部分を平坦化することができる。このような峰の平坦化処理は，平滑化処理で大きなフィルタを使用することでもできるが，その場合は，同時に輝度値変化領域３０の形状も平滑化され，樹冠の境界部分の輝度値変化が大きく損なわれるので，樹冠形状の検出をより困難にしてしまい好ましくない。
【００４１】
図１１は，上記膨張処理を１回行った後の画像であり，図１２は膨張処理を５回行った後の画像であり，更に，図１３は膨張処理を３４回行った後の画像である。膨張処理回数が多くなるほど，輝度の高い部分（明るい部分）の濃淡変化がなくなり平坦になっているのが理解される。
【００４２】
図１４は，縮退処理のフローチャート図であり，図１５は，縮退処理を説明する輝度値の空間変化のグラフ図である。また，図１６〜図１８は，膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。縮退処理は，平滑化処理された画像データの輝度値の空間変化について，その谷の部分を平坦化する処理である。図１５を参照しながら図１４のフローチャートを説明すると，まず，輝度値の空間変化の谷の深さを測定し，比較的多数の谷の深さを含む深さを所定値に設定する（Ｓ４０）。この所定値の考え方は，膨張処理と同じであり，樹冠の境界部分３０の輝度値の変化を変更することなく，適切に谷の部分を平坦化することができる値である。図１５の膨張処理後の輝度値の空間変化のグラフでは，谷部分は１つしかないので，その谷部分を平坦化できる深さが所定値として選択される。
【００４３】
次に，画素の輝度値を所定数加算した画像データを生成する（Ｓ４２）。つまり，図１５の実線の矢印のように，実線のグラフを破線のグラフのようにする。そして，元の輝度値を超えない範囲で近傍８画素の最小輝度値に置き換える処理を行う（Ｓ４４，Ｓ４６，Ｓ４８，Ｓ５０）。つまり，注目画素の輝度値が近傍８画素の最小輝度値よりも大きい場合は（Ｓ４４のＹＥＳ），更に、近傍８画素の最小輝度値が元の輝度値より小さくないなら（Ｓ４６のＮＯ）注目画素の輝度値を近傍８画素の最小輝度値に置き換えて（Ｓ４８），小さいなら（Ｓ４６のＹＥＳ）注目画素の輝度値を元の輝度値に置き換える（Ｓ５０）。そして，これらの置き換え処理Ｓ４４〜Ｓ５０を全ての画素について行う（Ｓ５２）。この置き換え処理は，複数回繰り返すことが望ましい。これが図１５の破線矢印の処理に該当する。破線矢印のように処理することで，一点鎖線に示した輝度値の画像データが生成される。
【００４４】
図１６は，上記縮退処理を１回行った後の画像であり，図１２は縮退処理を５回行った後の画像であり，更に，図１３は縮退処理を３４回行った後の画像である。縮退処理回数が多くなるほど，輝度の低い部分（暗い部分）の濃淡変化がなくなり平坦になっているのが理解される。
【００４５】
以上の膨張処理と縮退処理により平坦化された空間変化が求められる。上空から撮影した画像は，樹冠領域は光が反射し，輝度値が高くなるが，樹冠の間の領域は輝度値が低くなる傾向にある。但し，図１０に実線で示した平滑化後のグラフを利用して，前述の特許文献１に記載されたウォーターシェッド方法を適用して峰の部分を全て異なる樹冠として検出しようとすると，一つの樹冠内に複数の峰が存在する場合などはそれぞれが個別の領域と判定され樹冠が細分化される。そこで，上記のように平坦化処理を行うことで，本来なら一つである樹冠が細分化されてそれぞれが別の領域と判定されることを防ぐことができる。
【００４６】
図３に戻り，前述のとおり，膨張処理Ｓ１４Ａと縮退処理Ｓ１４Ｂとにより，画像データは樹冠の境界部分の輝度値変化領域３０を実質的に変化させることなく峰の部分と谷の部分とが平坦化される。次に，輝度値の空間変化に対する微分値が算出される（Ｓ１６Ａ）。
【００４７】
図１９は，領域分割処理Ｓ１６の微分値を示す図である。平坦化処理された画像データについて微分値を求めると，図１９に示されるとおり，平坦部分３２では，微分値はゼロになり（領域３８），輝度値が大きく変化する部分３０の微分値は大きくなる。なお，図１９には微分値の絶対値が示されている。また，輝度値が大きく変化する領域３０でもその傾きが一部変化する部分３４が存在する。輝度値の変化に段差が生じる部分である。この部分３４は，例えば，樹冠内において輝度値が変化する領域若しくは樹冠外において輝度値が変化する領域に対応する。そして，この傾きが変化する部分３４に対応する微分値は，極小値３６をとる。
【００４８】
前述のとおり，上空から撮影した画像は，樹冠の部分は光を反射しているので高い輝度値となり，樹冠以外の部分は影になり低い輝度値となる。したがって，画像の輝度値が峰の部分は樹冠の部分に対応し，谷の部分は樹冠以外の部分に対応する。そこで，図１９の輝度値のグラフにおいて，峰の部分の領域を樹冠領域とし，谷の部分を樹冠外領域と判定すればよいことになる。但し，微分値が極小値３６の部分３４が，樹冠内か樹冠外かを判定して，樹冠形状を正確に求めることが必要になる。
【００４９】
図２０は，微分値画像を示す図であり，微分値大を黒に微分値小を白で示した画像（１）と，微分値の極小箇所を黒にそれ以外を白で示した画像（２）とを示す。つまり，画像（１）は，輝度値が大きく変動する箇所が黒く示され，画像（２）は微分値がゼロの領域３８及び極小値の部分３６が黒く示される。従って，画像（２）の黒い部分は図１９の輝度値が平坦化された領域３２と輝度値の傾きが変動する部分３４とに対応し，白い部分が輝度値変動部分３０に対応する。したがって，平坦化領域３２の間に樹冠の境界線が位置することになる。但し，微分値が極小値になる部分３４を樹冠側にするか否かを判定する必要がある。
【００５０】
図３に戻り，ウォーターシェッドアルゴリズムによる領域分割処理Ｓ１６Ｂは，微分値データに対して，微分値が０と微分値が極小値をとる箇所を開始点とする領域拡張処理と，極小値から拡張した領域の併合処理とを有する。
【００５１】
図２１は，領域拡張処理と併合処理とを説明する図である。領域拡張処理では，図２１（１）に示されるとおり，微分値が０の箇所と極小値の箇所とにそれぞれ異なるラベル番号を割り当てる。それと共に，それぞれの領域には，輝度値が峰か谷かの属性データと，微分値０か微分値極小値かの属性データも割り当てられる。図２１の例では，領域Ｒ０〜Ｒ３の黒い領域にはそれぞれラベル番号０〜３が割り当てられ，属性データとして，領域Ｒ０とＲ１は微分値０，領域Ｒ２とＲ３は極小値が割り当てられ，領域Ｒ１には輝度値の峰，領域Ｒ０には輝度値の谷が割り当てられている。そして，ラベル番号が与えられた画素に隣接する画素にも同じラベル番号を割り当てるラベリング処理を繰り返す。このラベリング処理では，微分値の峰を越えては拡張することはできない。したがって，ラベリング処理による領域拡張を行うと，微分値の尾根が領域の境界となる。全ての画素にラベル番号が与えられるとこのラベリング処理は終了する。その結果，同じラベル番号がそれぞれ拡張領域Ｒ０〜Ｒ３となる。
【００５２】
図２１（１）に示すとおり，拡張処理後の領域には，微分値０から拡張した輝度値が峰の領域Ｒ１と，微分値０から拡張した輝度値が谷の領域Ｒ０と，微分値の極小値から拡張した領域Ｒ２およびＲ３とで構成される。図２２は，上記のウォーターシェッドアルゴリズムによる領域分割処理をした画像を示す図である。図２２（１）では，領域Ｒ１に対応する領域がグレーで示され，それ以外の領域Ｒ０，Ｒ２，Ｒ３に対応する領域は白で示されている。
【００５３】
次に，拡張領域のうち極小値から生成した領域Ｒ２とＲ３を輝度値の峰の領域Ｒ１か谷の領域Ｒ０かのいずれと同一領域なのかを判定し併合する処理が必要になる。そこで，極小値から生成した全ての領域Ｒ２とＲ３に対して，その領域の平均輝度値が最も近い隣接領域と同一領域であると判定し、その隣接領域に併合する。この併合処理では，極小値からの領域が輝度値の峰の領域Ｒ１あるいは谷の領域Ｒ０に併合される場合もあれば，隣接する極小値からの領域に併合される場合もある。したがって，上記の隣接領域への併合処理は，極小値からの領域がなくなるまで繰り返される。
【００５４】
図２１（２）には，領域Ｒ３が領域Ｒ１に併合され，領域Ｒ２が領域Ｒ０に併合された例が示されている。同様に，図２２（２）には，極小値からの領域が全て併合された画像が示されている。矢印５０が輝度値の峰の領域に併合され，矢印５２が輝度値の谷の領域に併合された例である。このように，極小値から領域拡張をして，平均輝度値の類似する領域に併合させると以下のメリットがある。図１９の平坦化された輝度値のグラフにおいて，単純に微分値０の領域を出発点として領域拡張を行ったのでは，微分値０の領域の中間位置が樹冠の境界になる。しかしながら，樹冠内側や外側の画像には輝度値が変化する部分があり，その部分を正しく把握して樹冠内部か外部かを判定しなければ，樹冠の形状を正確に検出することはできない。そこで，本実施の形態では，前述したとおり，微分値の極小値の点３６を出発点とする拡張領域を一旦生成し，その後，その平均輝度値に応じて隣接する領域に併合させている。これにより，より正確な樹冠形状を検出することができる。結局，図２１（２）に示される斜線の領域が樹冠領域に該当することになる。
【００５５】
［樹種の判定］
次に，樹冠の場所と形状が検出された後，各樹冠の樹種の判定が行われる。図２の評価方法に示すとおり，本実施の形態では，樹冠の画像（第１主成分の画像）の空間的模様（テクスチャ特徴量）を抽出し，樹種が既知の樹冠画像から求めた空間的模様（テクスチャ特徴量）からなるトレーニングデータと比較し，樹種を判定する（Ｓ１８）。ここで，テクスチャ特徴量としては，同時生起行列から求められるテクスチャ特徴量などである。そして，樹冠の位置と形状（大きさ）及び樹種を出力する（Ｓ２０）。そこで，樹種判定方法について詳述する。
【００５６】
図２３は，本実施の形態における樹冠の樹種判定方法のフローチャート図である。図２２（２）で求められた樹冠形状に基づき，図７の撮影データの第１主成分の画像について，同時生起行列を求める（Ｓ６０）。そして，同時生起行列からテクスチャ特徴量として局所一様性やコントラストなどを求める（Ｓ６２）。更に，実地調査などにより樹種が既知の樹冠画像についても同時生起行列を求め（Ｓ６４），その同時生起行列からテクスチャ特徴量を求め，それをトレーニングデータとする（Ｓ６６）。最後に，各樹冠のテクスチャ特徴量をトレーニングデータと比較し，一致するトレーニングデータの樹種に分類する（Ｓ６８）。
【００５７】
図２４は，同時生起行列を説明する図である。同時生起行列とは，テクスチャ特徴量のうち統計的なテクスチャ特徴の計算方法の一つであり，以下に定義される同時生起行列から画像の空間的模様の一つであるコントラストや局所一様性を求めることができる。図２４（１）に示される４×４画像を例にして説明すると，図２４（２）に示されるとおり，画像の輝度値ｉの画素から角度θの方向に長さｒの変位δ＝（ｒ，θ）離れた画素の輝度値がｊである確率Ｐδ（ｉ，ｊ）を要素とする行列が同時生起行列である。図２４（３）に距離ｒ＝１で角度θ＝０，４５，９０，１３５の４つの変位に対する同時生起行列を示す。変位（１，０°）の場合の同時生起行列を例にすると，横方向の隣（ｒ＝１）の画素間で輝度ｉ＝０，ｊ＝０になる個数は４回である。また，輝度ｉ＝０，ｊ＝１になる個数は２回であり，ｉ＝０，ｊ＝２になる個数は１回，ｉ＝０，ｊ＝３になる個数は０回である。同様に，ｉ＝１，ｊ＝０は２回，ｉ＝１，ｊ＝１は４回，ｉ＝１，ｊ＝２は０回，ｉ＝１，ｊ＝３も０回である。以下同様に横方向に隣接する画素間で輝度ｉ，ｊが所定の関係になる個数をそれぞれカウントすることで，変位（１，０°）の同時生起行列を求めることができる。更に，同時生起行列は，変位（１，４５°），（１，９０°），（１，１３５°）についても求められる。但し，図２４に示された要素の値とは異なり，要素の総和が１．０になるように正規化しておくことが望ましい。つまり，各要素が確率になる。
【００５８】
このように，同時生起行列は，ある画像（図２４（１））の画素間の輝度値の空間的配置関係を要素とする行列である。そして，この行列に対して，コントラストは，次の通りである。
【００５９】
コントラスト＝Σｋ2Ｐx-y（ｋ） （Σはｋ＝０〜（ｎ−１））
但し，Ｐx-y（ｋ）＝ΣΣＰδ（ｉ，ｊ） （│ｉ−ｊ│＝ｋ，ｋ＝０〜(n-1)，Σはｉ＝０〜(n-1)，ｊ＝０〜(n-1)）とする。
【００６０】
このように，コントラストは，│ｉ−ｊ│＝ｋの確率Ｐδ（ｉ，ｊ）を累積した確率Ｐx-y（ｋ）にｋ2を乗算した累積値であるので，同時生起行列の主対角から離れた（ｋが大きい画素対）確率がどのくらい高いかを示すものである。つまり，輝度差が大きい画素対が多く存在することを意味し，変位の角度θに直交する方向にエッジが存在する場合に大きな値になる。
【００６１】
局所一様性は，次の通りである。
【００６２】
局所一様性＝ΣΣ［１／｛１＋（ｉ−ｊ）2｝］Ｐδ（ｉ，ｊ） （Σはｉ＝０〜(n-1)，ｊ＝０〜(n-1)）
このように，局所一様性（inverse difference moment）は，確率Ｐδ（ｉ，ｊ）に対して主対角でより高くなる重み付けをした累積値であり，行列の主対角付近に高い確率が集中しているほど，つまり隣接画素が同じ輝度値をとる回数が多いほど高い値になるので，局所的に輝度値が一様になっている場合に高い値になる。したがって，局所一様性が高いほどコントラストは低くなり，局所一様性が低いほどコントラストは高くなる傾向にある。したがって，局所一様性をテクスチャ特徴量としてもよく，コントラストをテクスチャ特徴量にしてもよい。
【００６３】
図２４（１）の画像の場合，コントラストと局所一様性は，以下の値に計算される。
【００６４】
コントラスト：0.583 (0度), 1.000 (90度), 1.778 (45度), 0.444 (135度)
局所一様性：0.808 (0度), 0.700(90度), 0.511 (45度), 0.778 (135度)
統計的なテクスチャ特徴量の計算手法は，画像データに対してコンピュータソフトウエアによる情報処理を施すことにより空間的模様を自動的に求めることができ，大量の画像データに対する空間的特徴の抽出手法として適している。統計的なテクスチャ特徴量の計算方法として，上記の同時生起行列以外にも，（１）画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含むテクスチャ特徴量，（２）画像内で前記一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められるテクスチャ特徴量，（３）画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められるテクスチャ特徴量などでも，同様に空間的模様の指標として利用することができる。これらのテクスチャ特徴量については，前述の「新編，画像解析ハンドブック，高木幹雄・下田陽久監修」に詳しく説明されている。
【００６５】
本実施の形態では，上空から撮影した画像データは，単一画素あたり数１０センチメートル〜数メートルの分解能であり，単一樹冠内に複数行複数列の画素が含まれるので，樹冠の画像データから空間的模様を抽出することができる。したがって，単一の画素のスペクトル分析によるよりも空間的模様を基準にしたほうが樹冠の樹種をより適切に判定することができる。また，単一の画素内には複数の枝葉が含まれる程度の分解能であるので，樹冠の空間的模様は枝葉の模様とは異なる。
【００６６】
上記した本実施の形態での樹冠評価方法は，コンピュータソフトウエアによる情報処理により実現することができる。その場合は，汎用コンピュータに樹冠評価プログラムをインストールし，上空から撮影した画像データと樹種が既知の樹冠画像データとを入力し，情報処理することで，森林の樹冠の位置と形状を検出することができ，更に樹冠のテクスチャ特徴量を利用して樹種を判定することができる。
【図面の簡単な説明】
【００６７】
【図１】本実施の形態における人工衛星など上空から撮影した画像と樹冠との関係を示す図である。
【図２】本実施の形態における森林の樹冠の評価方法についてのフローチャート図である。
【図３】森林の樹冠形状の検出方法のフローチャート図である。
【図４】衛星画像の一例を示す図である。
【図５】衛星画像の一例を示す図である。
【図６】第１主成分の画像生成を説明するための図である。
【図７】図４，５の４つの画像データから求めた第１主成分の画像を示す図である。
【図８】図６の第１主成分の画像をガウシアンフィルタにより平滑化処理した画像を示す図である。
【図９】膨張処理のフローチャート図である。
【図１０】膨張処理を説明する輝度値の空間変化のグラフ図である。
【図１１】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１２】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１３】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１４】縮退処理のフローチャート図である。
【図１５】縮退処理を説明する輝度値の空間変化のグラフ図である。
【図１６】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１７】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１８】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１９】領域分割処理Ｓ１６の微分値を示す図である。
【図２０】微分値画像を示す図である。
【図２１】領域拡張処理と併合処理とを説明する図である。
【図２２】ウォーターシェッドアルゴリズムによる領域分割処理をした画像を示す図である。
【図２３】本実施の形態における樹冠の樹種判定方法のフローチャート図である。
【図２４】同時生起行列を説明する図である。
【符号の説明】
【００６８】
Ｓ１０〜Ｓ１６：樹冠検出工程 Ｓ１８，Ｓ２０：樹冠の樹種判定工程
【技術分野】
【０００１】
本発明は，森林の樹冠の評価方法及びその樹冠評価プログラムに関し，特に，衛星写真画像など森林を上空から撮影した高分解能の画像データに対して所定の画像処理を行って森林の樹冠の位置，形状（大きさ）を検出し樹種を判定する樹冠評価方法及びその樹冠評価プログラムに関する。
【背景技術】
【０００２】
森林地帯の衛星写真画像データに画像処理を加えることで森林の樹冠を評価することが注目されている。樹冠とは，森林を構成する樹木を上部から観察した時に，樹木の枝葉で構成された丸みを持った形状をいい，樹種とは樹木の種類をいう。森林を構成する樹冠の位置，形状，大きさの検出に加えてその樹木の種類（樹種）を判定することで，例えば二酸化炭素の排出権の評価や，森林の評価を行うことができる。
【０００３】
従来から，人工衛星や航空機から撮影した森林のカラー写真画像を画像処理して，森林の樹冠形状画像を求め，その樹冠形状画像と，実際に調査した樹冠に対応するカラー画像とを比較して，樹冠の樹種を判定する森林評価方法が提案されている。例えば，以下の特許文献１に記載されている。
【０００４】
この森林評価方法によれば，人工衛星や航空機で撮影した森林のカラーデジタル写真画像を白黒デジタル画像に変換し，平滑化処理し，樹冠の境界線の強調処理をし，ウォーターシェッドアルゴリズムにより樹冠形状を求め，樹冠形状画像内の画素のスペクトル分析により樹種を判定する。
【０００５】
また，人工衛星から撮影した森林地帯の写真画像から，樹木の活性度の指標である正規化植生指数（ＮＤＶＩ：Normalized Difference Vegetation Index）を算出する森林の評価方法や，森林の写真画像の濃淡パターンに対して，その一様性，方向性，コントラストなどを濃度共起行列を用いて求め，その特徴量から森林の樹種を判定する方法や，ヘリコプタから撮影した森林の写真画像に対して樹冠形状を近似円と仮定してその近似円を森林の写真画像に適合させて樹冠を検出する方法なども提案されている。例えば，以下の非特許文献１に記載されるとおりである。
【０００６】
更に，高分解能衛星，例えば地球観測衛星ＩＫＯＮＯＳ，の画像を利用した植生解析が日本スペースイメージング社から提案されている。これによれば，高分解能衛星画像を前処理し，領域分割して樹冠形状を求め，各分割領域のスペクトル解析を行って，その波長スペクトルを基にして樹種を分類することが提案されている。但し，領域分割の方法については公表されておらず，スペクトル解析の具体的方法も公表されていない。
【０００７】
そして，衛星ＩＫＯＮＯＳからの高分解能の森林の画像を現地調査で同定させた樹冠形状に重ねて，その樹冠の画素の平均値と標準偏差を計算し，その樹冠のスペクトル分析により樹種を分類することが提案されている。例えば，以下の非特許文献２に記載されるとおりである。
【特許文献１】特開２００１−３５７３８０号公報
【非特許文献１】「衛星リモートセンシングによる植生モニタリング−有峰地域の森林活性度調査−平成１２−１４年度共同研究成果報告書」平成１５年３月，金沢大学工学部村本研究室，北陸電力株式会社立地環境部
【非特許文献２】「高分解能ＩＫＯＮＯＳ衛星による針広混交林の樹種分類」加藤正人，森林航測 大１９８号 ６〜９ ２００２．１１
【発明の開示】
【発明が解決しようとする課題】
【０００８】
非特許文献１に記載されている衛星画像（例えばＬａｎｄｓａｔ）のような画素当たり３０ｍと低分解能の場合は，画素内に複数の樹冠が含まれ，そもそも樹冠形状を検出することができない。したがって，森林全体の評価をするに止まる。また，同文献に記載されているような航空機写真画像の場合は，樹冠検出は近似円による方法が提案されているにすぎず，また写真の濃淡パターンによる樹種判定では，樹木が密集した森林全体の樹種判定をするに止まっている。
【０００９】
更に，特許文献１に記載された衛星画像に平滑化処理等を行ってウォーターシェッドアルゴリズムにより樹冠形状を求める方法によると，輝度変化が最も小さい部分を中心にして，輝度の変化に沿って四方に領域を拡大させ，隣接する拡大領域が衝突する箇所を樹冠の境界線と定めるものであるが，樹冠の境界部分の輝度値変化を残す程度の平滑化処理だけでは、高分解能の衛星画像データの樹冠の輝度は単純な峰の形状にならないので，樹冠を正確に検出できない。また，樹種の判定に樹冠の画像のスペクトル解析を利用する方法が記載されているが，樹冠画像が複数の画素からなる高分解能の画像データの場合に波長スペクトルの相関性を利用した方法を適用すると，各画素が同じ波長スペクトルを有するとは限らず，正確に樹種を判定することは困難である。
【００１０】
また，非特許文献２では，高分解能の衛星写真画像を利用しているものの，樹冠形状の検出は，実地調査画像との重ね合わせによるものであり，実地調査自体が高コストであり現実的な方法とはいえない。また，スペクトル解析による樹種判定は，前述のとおり正確ではない。
【００１１】
そして，日本スペースイメージング社の提案では，前述のとおり樹冠形状の検出方法やスペクトル解析による樹種判定の具体的方法については記載されていない。
【００１２】
そこで，本発明の目的は，森林の高分解能の衛星写真画像から樹冠形状を正確に検出することができ，その樹冠の樹種を正確に判定することができる森林の樹冠評価方法及びその樹冠評価プログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記の目的を達成するために，本発明の第１の側面によれば，森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠評価方法において，前記画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行い，その平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，樹冠の形状を求める。更に，樹冠の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する。
【００１４】
本発明の第１の側面において，好ましい実施例では，平坦化処理が，前記画像データの輝度値を所定値だけ減算し元の輝度値を超えない範囲で近傍画素の最大輝度値に置き換える膨張処理と，画像データの輝度値を所定値だけ加算し元の輝度値を下回らない範囲で近傍画素の最小輝度値に置き換える縮退処理とを有する。この実施例によれば，コンピュータソフトウエアを利用した画像処理により膨張処理と縮退処理を簡単に行うことができ，それにより峰と谷が平坦化された画像データについて領域分割処理を簡単に行うことができる。
【００１５】
更に，好ましい実施例では，領域分割処理は，前記平坦化処理された画像データの微分値を求める微分値生成処理と，微分値ゼロや極小値の位置を中心に隣接する画素領域を拡張する領域拡張処理と，極小値の位置から拡張された領域を微分値ゼロの位置から拡張された領域に併合する併合処理とを有する。微分値を利用することにより，コンピュータソフトウエアを利用した画像処理で領域分割処理を簡単に行うことができる。
【００１６】
本発明の第１の側面において，好ましい実施例によれば，上空から撮影した画像データは，単一画素当たり数１０センチメートル〜数メートル程度の分解能であり，単一樹冠内に複数行複数列の画素が含まれ且つ単一画素内に複数の枝葉が含まれる程度の分解能である。かかる分解能を有する画像データの場合，樹冠画像には樹種に応じた固有の模様（テクスチャ）が生成されるので，スペクトル分析でなくテクスチャ特徴量による樹種判定を容易に行うことができる。
【００１７】
更に，好ましい実施例では，前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点の輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・・，ｎ−１）を要素とする同時生起行列から求められる第１のテクスチャ特徴量，画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含む第２のテクスチャ特徴量，画像内で一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められる第３のテクスチャ特徴量，画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められる第４のテクスチャ特徴量のいずれかを含む。これらのテクスチャ特徴量を利用することにより，適切に樹種を判定することができる。
【００１８】
上記の目的を達成するために，本発明の第２の側面は，第１の側面の評価方法をコンピュータに実行させる樹冠評価プログラムである。
【発明の効果】
【００１９】
上記本発明の第１の側面によれば，樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理を行うことにより，その後の領域分割処理において，ウォーターシェッドアルゴリズム（流域検出方法）による簡単な画像処理により樹冠形状の検出を正確に行うことができる。また，樹冠の画像のテクスチャ特徴量により樹種を判定するので，樹冠内の複数の画素のスペクトル分析よりもより正確に樹種を判定することができる。
【発明を実施するための最良の形態】
【００２０】
以下，図面にしたがって本発明の実施の形態について説明する。但し，本発明の技術的範囲はこれらの実施の形態に限定されず，特許請求の範囲に記載された事項とその均等物まで及ぶものである。
【００２１】
図１は，本実施の形態における人工衛星など上空から撮影した画像と樹冠との関係を示す図である。この図により本実施の形態が対象とする画像データの分解能を説明する。森林を上空から撮影すると，樹木の枝葉で構成された丸みを持った樹冠１０，１２，１４，１６が観察される。衛星ＩＫＯＮＯＳなどの高分解能衛星の画像データの場合は，その分解能が１画素当たり１メートル，または６０センチメートルであり，樹冠の直径が数メートルから数十メートルとすると，図示されるとおり各樹冠は複数行複数列の画素を含むことになる。図１の例では，実線の四角形が樹冠の画素に対応し，破線の四角形が樹冠以外の画素に対応する。また，１メートル／画素程度の分解能では，１画素内に樹木の複数の枝葉が含まれる。したがって，樹木の枝葉を複数の画素で表現するほどの高い分解能ではない。
【００２２】
このように，衛星ＩＫＯＮＯＳによる画像データは，樹冠の形状を目視判読できる程度に高い分解能を持つので，かかる画像データをコンピュータを利用して画像処理することによりその樹冠形状を検出することができる。そして，検出された樹冠形状の樹種を自動判定できれば，森林の評価方法としてその価値を高めることができる。近年，地球温暖化に伴う二酸化炭素の排出権ビジネスの始まりに際して，植林により創出される排出権を効率的且つ正確に評価するツールが求められているところであり，かかる衛星画像データを利用した自動森林評価方法は，排出権取引に欠かせないものである。特に，衛星画像データから樹木１本１本の形状を検出し，その樹種を判定できれば，正規化植生指数（ＮＤＶＩ：Normalized Difference Vegetation Index）による植生活性度（樹木の葉緑素の度合い）などの算出と組み合わせるなどして，より正確に二酸化炭素の吸収量を評価することができる。
【００２３】
図２は，本実施の形態における森林の樹冠の評価方法についてのフローチャート図である。この評価方法は，森林を上空から撮影した画像データをコンピュータプログラムにより画像処理を行うことで自動的に行うことができる。したがって，この評価方法は，汎用コンピュータに当該画像処理プログラムをインストールした評価コンピュータにより行うことができる。
【００２４】
また，図３は，上記評価方法における森林の樹冠形状の検出方法のフローチャート図である。これらのフローチャートにしたがい，且つ具体的画像例を参照しながら，森林の樹冠の評価方法を説明する。
【００２５】
図２には，森林樹冠の評価方法の概略が示されている。まず，衛星画像データを入力する（Ｓ１０）。この衛星画像データは，青，緑，赤，近赤外（Ｂ，Ｇ，Ｒ，ＮＩＲ）の各波長域の地表反射データからなるマルチスペクトル画像データである。図４，５に，これらの画像の一例を示す。図４には，青と緑の衛星画像が示される。また，図５には，赤と近赤外の衛星画像が示される。各画像データは，画素に対応する輝度値であり，例えば８ビット，２５６階調の輝度値データである。これらの衛星画像から目視判読されるとおり，複数の画素を含む樹冠が検出可能な程度の分解能を有している。
【００２６】
次に，４つの画像データについて主成分分析して，第１主成分の画像を生成する（Ｓ１２）。この工程については，後に詳述するが，要すれば，図４，５に示した４つの画像データを４次元の輝度値と考えて，それら４次元の輝度値から最も特徴の現れる軸（特徴の分散が最大になる軸）を求め，４つの画像データを線形変換して当該軸上の輝度値（第１主成分）を求めるものである。
【００２７】
そして，第１主成分の画像データに対して膨張処理と縮退処理を行い，画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更することなく，峰と谷の部分を削除して平坦にする（Ｓ１４）。この平坦化処理を行う前の第1主成分の画像データの輝度値の空間変化では，樹冠の境界部分の輝度値が大きく増減すると共に，樹冠内においても輝度値が増減する峰と谷が複数存在する。このような峰や谷は，ウォーターシェッドアルゴリズムによる領域分割処理を行うとそれぞれ個別の領域となり、一つの樹冠が細分化される。この細分化を防ぎ正確な樹冠形状を検出するためには，平坦化処理により峰や谷を平坦化する処理が有効である。
【００２８】
平坦化処理された画像データの輝度値の空間変化に対して，ウォーターシェッドアルゴリズム（流域検出方法）を適用して，樹冠の領域分割処理を行い，樹冠の位置や形状（領域）を求める（Ｓ１６）。ウォーターシェッドアルゴリズム自体は公知の手法であるが，そのうち，本実施の形態では，画像データの輝度値の空間変化について微分値を求め，微分値ゼロや微分値の極小値の位置を開始点として，周囲の隣接する画素領域を拡張する領域拡張処理を行う。平坦化された画像データの微分値を利用することで，ウォーターシェッドアルゴリズムによる領域分割処理をコンピュータプログラムにより容易に行うことができる。
【００２９】
樹冠形状が検出された後，その樹冠の画像データ，ここでは第１主成分の画像データ，の空間的模様であるテクスチャ特徴量を抽出し，既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較し，樹種を判定する（Ｓ１８）。本実施の形態では，衛星ＩＫＯＮＯＳの画像の分解能では樹冠の画像データが複数の画素により構成されることを考慮して，従来のスペクトル分析ではなく，テクスチャ特徴量を利用して樹種判定を行う。このテクスチャ特徴量は，例えば，同時生起行列から求められる局所一様性などの特徴量が有効である。
【００３０】
そして，最後に評価対象の森林における樹冠の位置，形状（大きさ），樹種を出力する（Ｓ２０）。この出力データと前述のＮＤＶＩ値を利用して，森林の排出権の評価を行うことができる。
【００３１】
［樹冠の検出］
次に，図３のフローチャートにしたがって，森林の樹冠の検出方法について詳述する。樹冠検出は，図２の工程Ｓ１０〜Ｓ１６により行われ，図３にはそれらの構成を詳述しており，同じ工程には同じ引用番号を与えている。まず，衛星画像データの入力（Ｓ１０）は，図２と同じである。次に，第１主成分の画像生成工程Ｓ１２が行われる。
【００３２】
図６は，第１主成分の画像生成を説明するための図である。第１主成分の画像を生成する主成分分析については，例えば，「新編，画像解析ハンドブック，高木幹雄・下田陽久監修」に詳しく説明されており，それによると，主成分分析とは，多変量の測定値（４色の輝度値）が与えられたとき，線形変換により変量間の相関をなくし，より少ない変量により測定対象の特徴を記述しようとする変換である。図６により簡単に説明すると，仮に２種の変量1ｘ，2ｘが測定された時，二次元座標（1ｘ軸と2ｘ軸）上には，黒点の変量２２がプロットされる。そして，これら変量の分散が最大になる破線の軸が主成分軸２０になる。
【００３３】
本実施の形態の４色の画像データの輝度値に適用して説明すると，図６の式（１）は，式（３）の条件の下で変量Ｚの分散が最大になるように，式（２）の係数ａを定めた場合の，第１主成分である変量Ｚを示す。この変量Ｚは，４色の輝度値である変量1ｘ，2ｘ，3ｘ，4ｘを係数ａ1,a2,a3,a4との一次結合で表したものである。
【００３４】
この第１主成分Ｚは，４色の輝度値について共分散行列の固有値方程式を解き，その最大固有値に対応する固有ベクトルを第１主成分ベクトルａとし，式（１）のとおり４つの輝度値1ｘ，2ｘ，3ｘ，4ｘに第１主成分ベクトルａを乗算して求める。図７は，図４，５の４つの画像データから求めた第１主成分の画像である。各画素の輝度値は，その分散が最大になっており，樹冠形状や樹冠の模様（テクスチャ特徴量）を求めるのに適した値に変換されている。
【００３５】
図３に戻り，次に，第１主成分の画像に対して，平滑化処理を行う（Ｓ１３）。この平滑化処理は，従来から知られているガウシアンフィルタを利用した処理である。ガウシアンフィルタの一例が図８に示されている。この例では，ガウシアンフィルタＧＦは３×３のフィルタであり，平滑化処理では，画像データの各画素の輝度値とその周囲８近傍の画素の輝度値に対して，このガウシアンフィルタＧＦの係数をそれぞれ乗算し，その合計値が中心の画素の輝度値に置き換えられる。これにより，第１主成分の画像データの空間変化における高周波成分が除去され平滑化される。このガウシアンフィルタＧＦは，８近傍の画素の輝度値を対応する係数に応じて加算するものであり，そのサイズを３×３と小さくすることで，樹冠を検出するための境界部分の輝度値変化を残しつつ，平滑化することができる。
【００３６】
図８は，図６の第１主成分の画像を上記ガウシアンフィルタにより平滑化処理した画像を示す図である。画像データの空間変化が平滑化されて，高周波成分が除去されているが，樹冠を検出するための境界部分の輝度値の大きな変化（低周波成分）は残されたままである。
【００３７】
次に，平坦化処理Ｓ１４のために，膨張処理により輝度値の空間変化の峰の部分を除去して平坦化し（Ｓ１４Ａ），縮退処理により輝度値の空間変化の谷の部分を除去して平坦化する（Ｓ１４Ｂ）。膨張処理と縮退処理はいずれを先に行っても良い。
【００３８】
図９は，膨張処理のフローチャート図であり，図１０は，膨張処理を説明する輝度値の空間変化のグラフ図である。また，図１１〜図１３は，図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。膨張処理は，平滑化処理された画像データの輝度値の空間変化について，その峰の部分を平坦化する処理である。図１０を参照しながら図９のフローチャートを説明すると，まず，輝度値の空間変化の峰の高さを測定し，比較的多数の峰の高さを含む高さを所定値に設定する（Ｓ２０）。図１０の実線は，平滑化処理された画像データの輝度値の空間変化を示す。峰とは極大値をとる部分であり，この例では４つの峰が存在している。この場合，峰Ｍ１とＭ２が一つの樹冠，峰Ｍ３が樹冠以外，峰Ｍ４が別の樹冠と予測される。そこで，膨張処理により，峰Ｍ１とＭ２を平坦化して一つにし，峰Ｍ３はなくし，峰４は上部の一部分のみを平坦化する。そのために，極大値の位置にある峰の高さを近傍の谷からの高さとして検出し，それらのうち比較的多数の峰の高さを含む高さを上記所定値と設定する。ここでは，峰Ｍ１，Ｍ２，Ｍ３のうちもっと高い高さに設定される。
【００３９】
次に，各画素の輝度値を前記所定値だけ減算した画像データを生成する（Ｓ２２）。その結果，図１０中の破線のグラフ（輝度値の空間変化）が得られる。この減算した画像データの輝度値を，元の輝度値を超えない範囲で８近傍の画素の最大輝度値に置き換える処理を，全ての画素に対して行う（Ｓ２４，Ｓ２６，Ｓ２８，Ｓ３０）。つまり，注目画素の輝度値が近傍８画素の最大輝度値よりも小さい場合は（Ｓ２４のＹＥＳ），更に，最大輝度値が元の輝度値より大きくないときは（Ｓ２６のＮＯ）注目画素の輝度値を最大輝度値に置き換え（Ｓ２８），最大輝度値が元の輝度値より大きいときは（Ｓ２６のＹＥＳ）注目画素の輝度値を元の輝度値に置き換える（Ｓ３０）。上記の処理を全ての画素について行う（Ｓ３２）。この置き換え処理Ｓ２４〜Ｓ３２は所定回数繰り返すことが望ましい。
【００４０】
その結果，図１０中の一点鎖線のグラフの輝度値が得られる。つまり，実線矢印のように減算処理された破線のグラフが，破線矢印のように，近くの峰の輝度値まで（但し元の輝度値以下）上昇され，最終的に小さな峰Ｍ１とＭ２が平坦化されて一つになり，Ｍ３が平坦化され，大きな峰Ｍ４の上部一部分が平坦化される。つまり，樹冠の境界部分の急激な輝度値変化領域３０を実質的に変更せず，小さな峰や大きな峰の上部一部分を平坦化することができる。このような峰の平坦化処理は，平滑化処理で大きなフィルタを使用することでもできるが，その場合は，同時に輝度値変化領域３０の形状も平滑化され，樹冠の境界部分の輝度値変化が大きく損なわれるので，樹冠形状の検出をより困難にしてしまい好ましくない。
【００４１】
図１１は，上記膨張処理を１回行った後の画像であり，図１２は膨張処理を５回行った後の画像であり，更に，図１３は膨張処理を３４回行った後の画像である。膨張処理回数が多くなるほど，輝度の高い部分（明るい部分）の濃淡変化がなくなり平坦になっているのが理解される。
【００４２】
図１４は，縮退処理のフローチャート図であり，図１５は，縮退処理を説明する輝度値の空間変化のグラフ図である。また，図１６〜図１８は，膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。縮退処理は，平滑化処理された画像データの輝度値の空間変化について，その谷の部分を平坦化する処理である。図１５を参照しながら図１４のフローチャートを説明すると，まず，輝度値の空間変化の谷の深さを測定し，比較的多数の谷の深さを含む深さを所定値に設定する（Ｓ４０）。この所定値の考え方は，膨張処理と同じであり，樹冠の境界部分３０の輝度値の変化を変更することなく，適切に谷の部分を平坦化することができる値である。図１５の膨張処理後の輝度値の空間変化のグラフでは，谷部分は１つしかないので，その谷部分を平坦化できる深さが所定値として選択される。
【００４３】
次に，画素の輝度値を所定数加算した画像データを生成する（Ｓ４２）。つまり，図１５の実線の矢印のように，実線のグラフを破線のグラフのようにする。そして，元の輝度値を超えない範囲で近傍８画素の最小輝度値に置き換える処理を行う（Ｓ４４，Ｓ４６，Ｓ４８，Ｓ５０）。つまり，注目画素の輝度値が近傍８画素の最小輝度値よりも大きい場合は（Ｓ４４のＹＥＳ），更に、近傍８画素の最小輝度値が元の輝度値より小さくないなら（Ｓ４６のＮＯ）注目画素の輝度値を近傍８画素の最小輝度値に置き換えて（Ｓ４８），小さいなら（Ｓ４６のＹＥＳ）注目画素の輝度値を元の輝度値に置き換える（Ｓ５０）。そして，これらの置き換え処理Ｓ４４〜Ｓ５０を全ての画素について行う（Ｓ５２）。この置き換え処理は，複数回繰り返すことが望ましい。これが図１５の破線矢印の処理に該当する。破線矢印のように処理することで，一点鎖線に示した輝度値の画像データが生成される。
【００４４】
図１６は，上記縮退処理を１回行った後の画像であり，図１２は縮退処理を５回行った後の画像であり，更に，図１３は縮退処理を３４回行った後の画像である。縮退処理回数が多くなるほど，輝度の低い部分（暗い部分）の濃淡変化がなくなり平坦になっているのが理解される。
【００４５】
以上の膨張処理と縮退処理により平坦化された空間変化が求められる。上空から撮影した画像は，樹冠領域は光が反射し，輝度値が高くなるが，樹冠の間の領域は輝度値が低くなる傾向にある。但し，図１０に実線で示した平滑化後のグラフを利用して，前述の特許文献１に記載されたウォーターシェッド方法を適用して峰の部分を全て異なる樹冠として検出しようとすると，一つの樹冠内に複数の峰が存在する場合などはそれぞれが個別の領域と判定され樹冠が細分化される。そこで，上記のように平坦化処理を行うことで，本来なら一つである樹冠が細分化されてそれぞれが別の領域と判定されることを防ぐことができる。
【００４６】
図３に戻り，前述のとおり，膨張処理Ｓ１４Ａと縮退処理Ｓ１４Ｂとにより，画像データは樹冠の境界部分の輝度値変化領域３０を実質的に変化させることなく峰の部分と谷の部分とが平坦化される。次に，輝度値の空間変化に対する微分値が算出される（Ｓ１６Ａ）。
【００４７】
図１９は，領域分割処理Ｓ１６の微分値を示す図である。平坦化処理された画像データについて微分値を求めると，図１９に示されるとおり，平坦部分３２では，微分値はゼロになり（領域３８），輝度値が大きく変化する部分３０の微分値は大きくなる。なお，図１９には微分値の絶対値が示されている。また，輝度値が大きく変化する領域３０でもその傾きが一部変化する部分３４が存在する。輝度値の変化に段差が生じる部分である。この部分３４は，例えば，樹冠内において輝度値が変化する領域若しくは樹冠外において輝度値が変化する領域に対応する。そして，この傾きが変化する部分３４に対応する微分値は，極小値３６をとる。
【００４８】
前述のとおり，上空から撮影した画像は，樹冠の部分は光を反射しているので高い輝度値となり，樹冠以外の部分は影になり低い輝度値となる。したがって，画像の輝度値が峰の部分は樹冠の部分に対応し，谷の部分は樹冠以外の部分に対応する。そこで，図１９の輝度値のグラフにおいて，峰の部分の領域を樹冠領域とし，谷の部分を樹冠外領域と判定すればよいことになる。但し，微分値が極小値３６の部分３４が，樹冠内か樹冠外かを判定して，樹冠形状を正確に求めることが必要になる。
【００４９】
図２０は，微分値画像を示す図であり，微分値大を黒に微分値小を白で示した画像（１）と，微分値の極小箇所を黒にそれ以外を白で示した画像（２）とを示す。つまり，画像（１）は，輝度値が大きく変動する箇所が黒く示され，画像（２）は微分値がゼロの領域３８及び極小値の部分３６が黒く示される。従って，画像（２）の黒い部分は図１９の輝度値が平坦化された領域３２と輝度値の傾きが変動する部分３４とに対応し，白い部分が輝度値変動部分３０に対応する。したがって，平坦化領域３２の間に樹冠の境界線が位置することになる。但し，微分値が極小値になる部分３４を樹冠側にするか否かを判定する必要がある。
【００５０】
図３に戻り，ウォーターシェッドアルゴリズムによる領域分割処理Ｓ１６Ｂは，微分値データに対して，微分値が０と微分値が極小値をとる箇所を開始点とする領域拡張処理と，極小値から拡張した領域の併合処理とを有する。
【００５１】
図２１は，領域拡張処理と併合処理とを説明する図である。領域拡張処理では，図２１（１）に示されるとおり，微分値が０の箇所と極小値の箇所とにそれぞれ異なるラベル番号を割り当てる。それと共に，それぞれの領域には，輝度値が峰か谷かの属性データと，微分値０か微分値極小値かの属性データも割り当てられる。図２１の例では，領域Ｒ０〜Ｒ３の黒い領域にはそれぞれラベル番号０〜３が割り当てられ，属性データとして，領域Ｒ０とＲ１は微分値０，領域Ｒ２とＲ３は極小値が割り当てられ，領域Ｒ１には輝度値の峰，領域Ｒ０には輝度値の谷が割り当てられている。そして，ラベル番号が与えられた画素に隣接する画素にも同じラベル番号を割り当てるラベリング処理を繰り返す。このラベリング処理では，微分値の峰を越えては拡張することはできない。したがって，ラベリング処理による領域拡張を行うと，微分値の尾根が領域の境界となる。全ての画素にラベル番号が与えられるとこのラベリング処理は終了する。その結果，同じラベル番号がそれぞれ拡張領域Ｒ０〜Ｒ３となる。
【００５２】
図２１（１）に示すとおり，拡張処理後の領域には，微分値０から拡張した輝度値が峰の領域Ｒ１と，微分値０から拡張した輝度値が谷の領域Ｒ０と，微分値の極小値から拡張した領域Ｒ２およびＲ３とで構成される。図２２は，上記のウォーターシェッドアルゴリズムによる領域分割処理をした画像を示す図である。図２２（１）では，領域Ｒ１に対応する領域がグレーで示され，それ以外の領域Ｒ０，Ｒ２，Ｒ３に対応する領域は白で示されている。
【００５３】
次に，拡張領域のうち極小値から生成した領域Ｒ２とＲ３を輝度値の峰の領域Ｒ１か谷の領域Ｒ０かのいずれと同一領域なのかを判定し併合する処理が必要になる。そこで，極小値から生成した全ての領域Ｒ２とＲ３に対して，その領域の平均輝度値が最も近い隣接領域と同一領域であると判定し、その隣接領域に併合する。この併合処理では，極小値からの領域が輝度値の峰の領域Ｒ１あるいは谷の領域Ｒ０に併合される場合もあれば，隣接する極小値からの領域に併合される場合もある。したがって，上記の隣接領域への併合処理は，極小値からの領域がなくなるまで繰り返される。
【００５４】
図２１（２）には，領域Ｒ３が領域Ｒ１に併合され，領域Ｒ２が領域Ｒ０に併合された例が示されている。同様に，図２２（２）には，極小値からの領域が全て併合された画像が示されている。矢印５０が輝度値の峰の領域に併合され，矢印５２が輝度値の谷の領域に併合された例である。このように，極小値から領域拡張をして，平均輝度値の類似する領域に併合させると以下のメリットがある。図１９の平坦化された輝度値のグラフにおいて，単純に微分値０の領域を出発点として領域拡張を行ったのでは，微分値０の領域の中間位置が樹冠の境界になる。しかしながら，樹冠内側や外側の画像には輝度値が変化する部分があり，その部分を正しく把握して樹冠内部か外部かを判定しなければ，樹冠の形状を正確に検出することはできない。そこで，本実施の形態では，前述したとおり，微分値の極小値の点３６を出発点とする拡張領域を一旦生成し，その後，その平均輝度値に応じて隣接する領域に併合させている。これにより，より正確な樹冠形状を検出することができる。結局，図２１（２）に示される斜線の領域が樹冠領域に該当することになる。
【００５５】
［樹種の判定］
次に，樹冠の場所と形状が検出された後，各樹冠の樹種の判定が行われる。図２の評価方法に示すとおり，本実施の形態では，樹冠の画像（第１主成分の画像）の空間的模様（テクスチャ特徴量）を抽出し，樹種が既知の樹冠画像から求めた空間的模様（テクスチャ特徴量）からなるトレーニングデータと比較し，樹種を判定する（Ｓ１８）。ここで，テクスチャ特徴量としては，同時生起行列から求められるテクスチャ特徴量などである。そして，樹冠の位置と形状（大きさ）及び樹種を出力する（Ｓ２０）。そこで，樹種判定方法について詳述する。
【００５６】
図２３は，本実施の形態における樹冠の樹種判定方法のフローチャート図である。図２２（２）で求められた樹冠形状に基づき，図７の撮影データの第１主成分の画像について，同時生起行列を求める（Ｓ６０）。そして，同時生起行列からテクスチャ特徴量として局所一様性やコントラストなどを求める（Ｓ６２）。更に，実地調査などにより樹種が既知の樹冠画像についても同時生起行列を求め（Ｓ６４），その同時生起行列からテクスチャ特徴量を求め，それをトレーニングデータとする（Ｓ６６）。最後に，各樹冠のテクスチャ特徴量をトレーニングデータと比較し，一致するトレーニングデータの樹種に分類する（Ｓ６８）。
【００５７】
図２４は，同時生起行列を説明する図である。同時生起行列とは，テクスチャ特徴量のうち統計的なテクスチャ特徴の計算方法の一つであり，以下に定義される同時生起行列から画像の空間的模様の一つであるコントラストや局所一様性を求めることができる。図２４（１）に示される４×４画像を例にして説明すると，図２４（２）に示されるとおり，画像の輝度値ｉの画素から角度θの方向に長さｒの変位δ＝（ｒ，θ）離れた画素の輝度値がｊである確率Ｐδ（ｉ，ｊ）を要素とする行列が同時生起行列である。図２４（３）に距離ｒ＝１で角度θ＝０，４５，９０，１３５の４つの変位に対する同時生起行列を示す。変位（１，０°）の場合の同時生起行列を例にすると，横方向の隣（ｒ＝１）の画素間で輝度ｉ＝０，ｊ＝０になる個数は４回である。また，輝度ｉ＝０，ｊ＝１になる個数は２回であり，ｉ＝０，ｊ＝２になる個数は１回，ｉ＝０，ｊ＝３になる個数は０回である。同様に，ｉ＝１，ｊ＝０は２回，ｉ＝１，ｊ＝１は４回，ｉ＝１，ｊ＝２は０回，ｉ＝１，ｊ＝３も０回である。以下同様に横方向に隣接する画素間で輝度ｉ，ｊが所定の関係になる個数をそれぞれカウントすることで，変位（１，０°）の同時生起行列を求めることができる。更に，同時生起行列は，変位（１，４５°），（１，９０°），（１，１３５°）についても求められる。但し，図２４に示された要素の値とは異なり，要素の総和が１．０になるように正規化しておくことが望ましい。つまり，各要素が確率になる。
【００５８】
このように，同時生起行列は，ある画像（図２４（１））の画素間の輝度値の空間的配置関係を要素とする行列である。そして，この行列に対して，コントラストは，次の通りである。
【００５９】
コントラスト＝Σｋ2Ｐx-y（ｋ） （Σはｋ＝０〜（ｎ−１））
但し，Ｐx-y（ｋ）＝ΣΣＰδ（ｉ，ｊ） （│ｉ−ｊ│＝ｋ，ｋ＝０〜(n-1)，Σはｉ＝０〜(n-1)，ｊ＝０〜(n-1)）とする。
【００６０】
このように，コントラストは，│ｉ−ｊ│＝ｋの確率Ｐδ（ｉ，ｊ）を累積した確率Ｐx-y（ｋ）にｋ2を乗算した累積値であるので，同時生起行列の主対角から離れた（ｋが大きい画素対）確率がどのくらい高いかを示すものである。つまり，輝度差が大きい画素対が多く存在することを意味し，変位の角度θに直交する方向にエッジが存在する場合に大きな値になる。
【００６１】
局所一様性は，次の通りである。
【００６２】
局所一様性＝ΣΣ［１／｛１＋（ｉ−ｊ）2｝］Ｐδ（ｉ，ｊ） （Σはｉ＝０〜(n-1)，ｊ＝０〜(n-1)）
このように，局所一様性（inverse difference moment）は，確率Ｐδ（ｉ，ｊ）に対して主対角でより高くなる重み付けをした累積値であり，行列の主対角付近に高い確率が集中しているほど，つまり隣接画素が同じ輝度値をとる回数が多いほど高い値になるので，局所的に輝度値が一様になっている場合に高い値になる。したがって，局所一様性が高いほどコントラストは低くなり，局所一様性が低いほどコントラストは高くなる傾向にある。したがって，局所一様性をテクスチャ特徴量としてもよく，コントラストをテクスチャ特徴量にしてもよい。
【００６３】
図２４（１）の画像の場合，コントラストと局所一様性は，以下の値に計算される。
【００６４】
コントラスト：0.583 (0度), 1.000 (90度), 1.778 (45度), 0.444 (135度)
局所一様性：0.808 (0度), 0.700(90度), 0.511 (45度), 0.778 (135度)
統計的なテクスチャ特徴量の計算手法は，画像データに対してコンピュータソフトウエアによる情報処理を施すことにより空間的模様を自動的に求めることができ，大量の画像データに対する空間的特徴の抽出手法として適している。統計的なテクスチャ特徴量の計算方法として，上記の同時生起行列以外にも，（１）画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含むテクスチャ特徴量，（２）画像内で前記一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められるテクスチャ特徴量，（３）画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められるテクスチャ特徴量などでも，同様に空間的模様の指標として利用することができる。これらのテクスチャ特徴量については，前述の「新編，画像解析ハンドブック，高木幹雄・下田陽久監修」に詳しく説明されている。
【００６５】
本実施の形態では，上空から撮影した画像データは，単一画素あたり数１０センチメートル〜数メートルの分解能であり，単一樹冠内に複数行複数列の画素が含まれるので，樹冠の画像データから空間的模様を抽出することができる。したがって，単一の画素のスペクトル分析によるよりも空間的模様を基準にしたほうが樹冠の樹種をより適切に判定することができる。また，単一の画素内には複数の枝葉が含まれる程度の分解能であるので，樹冠の空間的模様は枝葉の模様とは異なる。
【００６６】
上記した本実施の形態での樹冠評価方法は，コンピュータソフトウエアによる情報処理により実現することができる。その場合は，汎用コンピュータに樹冠評価プログラムをインストールし，上空から撮影した画像データと樹種が既知の樹冠画像データとを入力し，情報処理することで，森林の樹冠の位置と形状を検出することができ，更に樹冠のテクスチャ特徴量を利用して樹種を判定することができる。
【図面の簡単な説明】
【００６７】
【図１】本実施の形態における人工衛星など上空から撮影した画像と樹冠との関係を示す図である。
【図２】本実施の形態における森林の樹冠の評価方法についてのフローチャート図である。
【図３】森林の樹冠形状の検出方法のフローチャート図である。
【図４】衛星画像の一例を示す図である。
【図５】衛星画像の一例を示す図である。
【図６】第１主成分の画像生成を説明するための図である。
【図７】図４，５の４つの画像データから求めた第１主成分の画像を示す図である。
【図８】図６の第１主成分の画像をガウシアンフィルタにより平滑化処理した画像を示す図である。
【図９】膨張処理のフローチャート図である。
【図１０】膨張処理を説明する輝度値の空間変化のグラフ図である。
【図１１】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１２】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１３】図８の平滑化処理された画像データに膨張処理を行った結果の画像例を示す図である。
【図１４】縮退処理のフローチャート図である。
【図１５】縮退処理を説明する輝度値の空間変化のグラフ図である。
【図１６】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１７】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１８】膨張処理された画像データに縮退処理を行った結果の画像例を示す図である。
【図１９】領域分割処理Ｓ１６の微分値を示す図である。
【図２０】微分値画像を示す図である。
【図２１】領域拡張処理と併合処理とを説明する図である。
【図２２】ウォーターシェッドアルゴリズムによる領域分割処理をした画像を示す図である。
【図２３】本実施の形態における樹冠の樹種判定方法のフローチャート図である。
【図２４】同時生起行列を説明する図である。
【符号の説明】
【００６８】
Ｓ１０〜Ｓ１６：樹冠検出工程 Ｓ１８，Ｓ２０：樹冠の樹種判定工程
【特許請求の範囲】
【請求項１】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価方法において，
前記撮影した画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理工程と，
当該平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，前記森林内の樹冠の形状を求める工程とを有することを特徴とする樹冠の評価方法。
【請求項２】
請求項１において，
前記平坦化処理工程は，前記画像データの輝度値を第１の所定値だけ減算し元の輝度値を超えない範囲で近傍画素の最大輝度値に置き換える膨張処理工程と，前記画像データの輝度値を第２の所定値だけ加算し元の輝度値を下回らない範囲で近傍画素の最小輝度値に置き換える縮退処理工程とを有することを特徴とする樹冠の評価方法。
【請求項３】
請求項１において，
前記領域分割処理は，前記平坦化処理された画像データの微分値を求める微分値生成処理工程と，微分値ゼロ及び極小値の位置を中心に隣接する画素領域を拡張する領域拡張処理工程と，前記極小値の位置から拡張された領域を前記微分値ゼロの位置から拡張された領域に併合する併合処理とを有することを特徴とする樹冠の評価方法。
【請求項４】
請求項１において，
更に，前記上空から撮影した画像データは，複数の色に対応する多次元の輝度値であり，当該多次元の輝度値の分散共分散行列から固有ベクトルを求め，当該多次元の輝度値に前記固有ベクトルを乗算して第１主成分の画像データを求める第１主成分生成工程を有し，
前記平坦化処理は，前記撮像した画像データに代えて当該第１主成分の画像データに対して行われることを特徴とする樹冠の評価方法。
【請求項５】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価方法において，
前記撮影した画像データを画像処理して前記森林内の樹冠の形状を求める工程と，
前記求められた樹冠形状の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する樹種判定工程とを有することを特徴とする樹冠の評価方法。
【請求項６】
請求項５において，
前記上空から撮影した画像データは，単一画素当たり数１０センチメートル〜数メートルの分解能であり，単一樹冠内に複数行複数列の画素が含まれ且つ単一画素内に複数の枝葉が含まれる程度の分解能であることを特徴とする樹冠の評価方法。
【請求項７】
請求項５または６において，
前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点が輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・，ｎ−１）を要素とする同時生起行列から求められる第１のテクスチャ特徴量，画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含む第２のテクスチャ特徴量，画像内で前記一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められる第３のテクスチャ特徴量，画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められる第４のテクスチャ特徴量のいずれかを含むことを特徴とする樹冠の評価方法。
【請求項８】
請求項６において，
前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点の輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・，ｎ−１）を要素とする同時生起行列を求め，当該同時生起行列から求められるコントラスト，局所一様性のいずれかを含むことを特徴とする樹冠の評価方法。
【請求項９】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価手順をコンピュータに実行させる樹冠評価プログラムにおいて，前記評価手順は，
前記撮影した画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理手順と，
当該平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，前記森林内の樹冠の形状を求める手順とを有することを特徴とする樹冠の評価プログラム。
【請求項１０】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価手順をコンピュータに実行させる樹冠評価プログラムにおいて，前記評価手順は，
前記撮影した画像データを画像処理して前記森林内の樹冠の形状を求める手順と，
前記求められた樹冠形状の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する樹種判定手順とを有することを特徴とする樹冠の評価プログラム。
【請求項１】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価方法において，
前記撮影した画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理工程と，
当該平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，前記森林内の樹冠の形状を求める工程とを有することを特徴とする樹冠の評価方法。
【請求項２】
請求項１において，
前記平坦化処理工程は，前記画像データの輝度値を第１の所定値だけ減算し元の輝度値を超えない範囲で近傍画素の最大輝度値に置き換える膨張処理工程と，前記画像データの輝度値を第２の所定値だけ加算し元の輝度値を下回らない範囲で近傍画素の最小輝度値に置き換える縮退処理工程とを有することを特徴とする樹冠の評価方法。
【請求項３】
請求項１において，
前記領域分割処理は，前記平坦化処理された画像データの微分値を求める微分値生成処理工程と，微分値ゼロ及び極小値の位置を中心に隣接する画素領域を拡張する領域拡張処理工程と，前記極小値の位置から拡張された領域を前記微分値ゼロの位置から拡張された領域に併合する併合処理とを有することを特徴とする樹冠の評価方法。
【請求項４】
請求項１において，
更に，前記上空から撮影した画像データは，複数の色に対応する多次元の輝度値であり，当該多次元の輝度値の分散共分散行列から固有ベクトルを求め，当該多次元の輝度値に前記固有ベクトルを乗算して第１主成分の画像データを求める第１主成分生成工程を有し，
前記平坦化処理は，前記撮像した画像データに代えて当該第１主成分の画像データに対して行われることを特徴とする樹冠の評価方法。
【請求項５】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価方法において，
前記撮影した画像データを画像処理して前記森林内の樹冠の形状を求める工程と，
前記求められた樹冠形状の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する樹種判定工程とを有することを特徴とする樹冠の評価方法。
【請求項６】
請求項５において，
前記上空から撮影した画像データは，単一画素当たり数１０センチメートル〜数メートルの分解能であり，単一樹冠内に複数行複数列の画素が含まれ且つ単一画素内に複数の枝葉が含まれる程度の分解能であることを特徴とする樹冠の評価方法。
【請求項７】
請求項５または６において，
前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点が輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・，ｎ−１）を要素とする同時生起行列から求められる第１のテクスチャ特徴量，画像の輝度値に対し全体が１．０になるよう正規化された輝度値ヒストグラムＰ（ｉ）（ｉは輝度値）の平均，分散，歪度，尖度を含む第２のテクスチャ特徴量，画像内で前記一定の変位δ＝（ｒ，θ）だけ離れた２点の輝度値の差がｋである確率Ｐδ（ｋ）（ｋ＝０，１，・・・ｎ−１）から求められる第３のテクスチャ特徴量，画像内でθ方向の輝度値ｉの点がｊ個続く頻度Ｐθ（ｉ，ｊ）を要素とするランレングス行列から求められる第４のテクスチャ特徴量のいずれかを含むことを特徴とする樹冠の評価方法。
【請求項８】
請求項６において，
前記テクスチャ特徴量は，画像の輝度値ｉの点から角度θの方向に長さｒの変位δ＝（ｒ，θ）だけ離れた点の輝度値ｊである確率Ｐδ（ｉ，ｊ）（ｉ，ｊ＝０，１，・・・，ｎ−１）を要素とする同時生起行列を求め，当該同時生起行列から求められるコントラスト，局所一様性のいずれかを含むことを特徴とする樹冠の評価方法。
【請求項９】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価手順をコンピュータに実行させる樹冠評価プログラムにおいて，前記評価手順は，
前記撮影した画像データの輝度値の空間変化について樹冠の境界部分の輝度値変化を実質的に変更せずに峰と谷の部分を平坦にする平坦化処理手順と，
当該平坦化された画像データの前記輝度値の空間変化に対して領域分割処理を行って，前記森林内の樹冠の形状を求める手順とを有することを特徴とする樹冠の評価プログラム。
【請求項１０】
森林を上空から撮影した画像データを処理して森林の樹冠を評価する樹冠の評価手順をコンピュータに実行させる樹冠評価プログラムにおいて，前記評価手順は，
前記撮影した画像データを画像処理して前記森林内の樹冠の形状を求める手順と，
前記求められた樹冠形状の画像データのテクスチャ特徴量を求め，当該樹冠の画像のテクスチャ特徴量を既知の樹種の樹冠画像データから求めた基準テクスチャ特徴量と比較して樹冠の樹種を判定する樹種判定手順とを有することを特徴とする樹冠の評価プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【公開番号】特開２００６−２８５３１０（Ｐ２００６−２８５３１０Ａ）
【公開日】平成１８年１０月１９日（２００６．１０．１９）
【国際特許分類】
【出願番号】特願２００５−１００４７０（Ｐ２００５−１００４７０）
【出願日】平成１７年３月３１日（２００５．３．３１）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１７年２月９日 金沢大学主催の「自然科学研究科電子情報システム専攻 平成１６年度 修士論文口頭発表会」において文書をもって発表
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１７年２月２２日 金沢大学主催の「平成１６年度 卒業研究発表会」において文書をもって発表
【出願人】（５０４１６０７８１）国立大学法人金沢大学 (282)
【出願人】（０００２４２６４４）北陸電力株式会社 (112)
【Ｆターム（参考）】
【公開日】平成１８年１０月１９日（２００６．１０．１９）
【国際特許分類】
【出願日】平成１７年３月３１日（２００５．３．３１）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１７年２月９日 金沢大学主催の「自然科学研究科電子情報システム専攻 平成１６年度 修士論文口頭発表会」において文書をもって発表
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１７年２月２２日 金沢大学主催の「平成１６年度 卒業研究発表会」において文書をもって発表
【出願人】（５０４１６０７８１）国立大学法人金沢大学 (282)
【出願人】（０００２４２６４４）北陸電力株式会社 (112)
【Ｆターム（参考）】
[ Back to top ]