JP2015023858A - 林相解析装置、林相解析方法及びプログラム - Google Patents

Abstract

【課題】航空写真等から森林の変化領域を高精度で自動的に生成することを可能とする。
【解決手段】統計分布作成部２６は、森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、各林相の注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める。変化領域検出部２８は各注目領域での統計分布における外れ値を検出し、当該外れ値を与える要素領域を当該注目領域において第１の時期との間で林相が変化した変化領域として検出する。
【選択図】図１

Description

本発明は、上空から森林を撮影した画像と当該画像とは別の時期の既存の林相区分図とから森林の変化領域を抽出する林相解析装置、林相解析方法及びプログラムに関する。

航空機などによる上空からのリモートセンシングは、広範囲での地上の状況把握が可能であり、伐採や風倒被害等による森林変化の検出や、森林の林相区分図の作成に利用されている。ちなみに、我が国では一般に数年を単位として森林の計画・管理が行われており、これに対応して森林の調査や林相区分図の更新も数年程度ごとに行われる。

既存の林相区分図と新しく取得した航空写真等の高分解能画像データとを利用した森林変化の検出や林相区分図の更新は従来、基本的に人間が目視判読によって行っている。しかしながら、広範囲にわたる目視判読には膨大な時間を要するため、より効率的な森林変化検出、並びに林相区分図更新の手法が求められている。

森林変化を自動判読する従来技術として、新年度の画像から作成された林相区分図を既存の林相区分図と比較する手法が知られている。また、森林を撮影した画像から林相区分図を自動生成する技術に関しては、当該画像を領域分割して、類似した特徴を持つ複数の小領域である林相区画を生成し、当該林相区画ごとに林相分類を行う手法が研究されている。

特開２００８−０４６８３７号公報 特開２０１１−０２４４７１号公報

森林変化の自動検出において、既存の林相区分図と比較される新年度の林相区分図を自動生成する際に、森林を撮影した画像から取得した教師データを用いて林相分類が行われる。従来、教師データの取得に際しては人間が各林相に該当する領域を選択し、当該領域から各林相の教師データが抽出される。そのため、林相区分図の作成は完全には自動化されておらず、森林変化の検出作業も十分に効率化・自動化されていなかった。

また、上述の林相に関する教師データの取得は人間の作業負荷に起因して画像の一部分から行われる。ここで、同じ林相であってもその特徴は画像内でばらつきを有する。そのため、画像の一部分から得られた教師データを用いて生成した林相区分図には誤分類が生じやすい。そのため、自動生成された林相区分図を直接利用して森林の変化検出を行うと誤検出が発生しやすいという問題があった。

本発明は上記問題点を解決するためになされたものであり、既存の林相区分図と航空写真等の高分解能の画像とを利用して森林の変化領域を高精度で自動的に検出することを可能とし、ひいては林相区分図の自動生成を可能とする林相解析装置、林相解析方法及びプログラムを提供することを目的とする。

（１）本発明に係る林相解析装置は、森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成手段と、前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出手段と、を有する。

（２）上記（１）に記載する林相解析装置において、さらに、前記第１の時期の前記林相区分図に前記変化領域を反映させて前記第２の時期の林相区分図を生成する林相区分図生成手段を有する構成とすることもできる。

（３）上記（１）又は（２）に記載する林相解析装置において、さらに、前記原画像の１又は複数種類の区画用特徴量の類似性に基づいて前記対象地域を複数領域に分割して前記森林の林相区画を生成し出力する林相区画生成手段を備え、前記統計分布作成手段及び前記変化領域検出手段における前記要素領域が前記林相区画である構成とすることもできる。

（４）上記（３）に記載する林相解析装置において、さらに、前記原画像から、当該原画像の隣接する複数画素からなる領域を画素として分解能を低下させた劣化画像を生成する画像劣化手段を備え、前記林相区画生成手段が、前記劣化画像を、その１又は複数種類の前記区画用特徴量の類似性に基づいて複数領域に分割し林相区画を生成する初期区画手段と、隣接する林相区画をそれらの１又は複数種類の前記区画用特徴量の類似性に基づいて結合して新たな林相区画を生成する逐次区画処理を１回又は複数回行い、前記対象地域における前記森林の林相区画を出力する後続区画手段と、を備え、前記後続区画手段は少なくとも１回の前記逐次区画処理を前記原画像の前記区画用特徴量、又は前記原画像及び前記劣化画像それぞれの前記区画用特徴量に基づいて行う構成とすることができる。

（５）上記（３）又は（４）に記載する林相解析装置において、前記原画像からテクスチャ情報を抽出するテクスチャ情報抽出手段を備え、前記統計分布作成手段が、前記林相区画生成手段から出力される前記各林相区画から抽出された前記テクスチャ情報を前記特徴情報として前記統計分布を求める構成とすることができる。

（６）上記（５）に記載する林相解析装置において、前記原画像からスペクトル情報を抽出するスペクトル情報抽出手段を備え、前記統計分布作成手段が、前記林相区画生成手段から出力される前記各林相区画から抽出された前記テクスチャ情報及び前記スペクトル情報を前記特徴情報として前記統計分布を求める構成とすることができる。

（７）上記（３）又は（４）に記載する林相解析装置において、前記原画像からスペクトル情報を抽出するスペクトル情報抽出手段を有し、前記統計分布作成手段は、前記林相区画生成手段から出力される前記各林相区画から抽出された前記スペクトル情報を前記特徴情報として前記統計分布を求める構成とすることができる。

（８）本発明に係る林相解析方法は、森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成ステップと、前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出ステップと、を有する。

（９）本発明に係るプログラムは、コンピュータに、森林を含む対象地域における林相解析を行わせるためのプログラムであって、当該コンピュータを、森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成手段、及び、前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出手段、として機能させる。

本発明によれば、既存の林相区分図と航空写真等の高分解能の画像とを利用して森林の変化領域を高精度で自動的に検出することが可能となり、ひいては林相区分図の自動生成が可能となる。

本発明の第１の実施形態である林相解析システムの概略の構成を示すブロック図である。 日向／日陰分布パターンの例を示す説明図である。 日向／日陰境界パターンの例を示す説明図である。 ＬＢＰ画像の例を示す説明図である。 正規化処理前のオルソ画像の例を示す画像である。 図５の画像を正規化処理して得られた画像である。 本発明の実施形態である林相解析システムにおける概略のデータフロー図である。 林相区画生成部による林相区画の生成処理を説明する概略のフロー図である。 階層的な領域分割処理の最初の階層にて低分解能画像から生成された林相区画の例を示す画像である。 階層的な領域分割処理の途中階層にて生成された林相区画の例を示す画像である。 図１０に示す階層より後の階層にて生成された林相区画の例を示す画像である。 テクスチャ解析部及びスペクトル解析部による特徴情報データの抽出処理を説明する概略のフロー図である。 注目領域に関する統計分布の模式図である。 本発明の第２の実施形態である林相解析システムの概略の構成を示すブロック図である。

以下、本発明の実施の形態（以下実施形態という）について、図面に基づいて説明する。

［第１の実施形態］
図１は、第１の実施形態である林相解析システム２の概略の構成を示すブロック図である。本システムは、演算処理装置４、記憶装置６、入力装置８及び出力装置１０を含んで構成される。演算処理装置４として、本システムの処理を行う専用のハードウェアを作ることも可能であるが、本実施形態では演算処理装置４は、コンピュータ及び、当該コンピュータ上で実行されるプログラムを用いて構築される。

演算処理装置４は、コンピュータのＣＰＵ（Central Processing Unit）からなり、後述する林相区画生成部２０、テクスチャ解析部２２、スペクトル解析部２４、統計分布作成部２６、変化領域検出部２８、林相区分図生成部３０として機能する。

記憶装置６はＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等のメモリ装置である。記憶装置６は演算処理装置４にて実行される各種のプログラムや、本システムの処理に必要な各種データを記憶し、演算処理装置４との間でこれらの情報を入出力する。例えば、記憶装置６には、オルソ画像データ４０、林相区分図データ４２、林相区画に関するスペクトル条件４４及び幾何条件４６が予め格納される。

オルソ画像データ４０は航空機等から撮影された空中写真画像に基づいて生成される。本システムではオルソ画像データ４０は、森林を含む対象地域を上空から撮影した原画像として用いられる高分解能画像データであり、例えば、赤（Ｒ）、緑（Ｇ）、青（Ｂ）の３成分からなる３バンドのマルチスペクトル画像、又はこれらに近赤外（ＮＩＲ）を加えた４成分からなる４バンドのマルチスペクトル画像である。

林相区分図データ４２は、オルソ画像データ４０に映る対象地域についての既存の林相区分図を表すデータである。林相区分図データ４２は対象地域を領域分割して得られる要素領域である林相区画と、当該区画における林相とを表す。ここで、林相区分図データ４２は時期Ｔ１に関するものであり、一方、オルソ画像データ４０は時期Ｔ１とは異なる時期Ｔ２に関するものであるとする。Ｔ１とＴ２との前後関係は任意であるが、本実施形態に示す例ではＴ１がＴ２に先行するものとする。

スペクトル条件４４及び幾何条件４６は、原画像を領域分割して生成されるオブジェクト領域である林相区画が満たすべき条件である。スペクトル条件４４は林相区画を構成する画素群の色（Color）や濃度・明るさから抽出されるスペクトル特徴量に関する条件を規定する。例えば、スペクトル特徴量として各色（バンド）の平均値、標準偏差などの統計量を用いることができる。一方、幾何条件４６は林相区画の幾何特徴量である形状（Shape）に関する条件を規定する。例えば、形状は、コンパクトネス（Compactness）やスムースネス（Smoothness）といったパラメータで表される。

入力装置８は、キーボード、マウスなどであり、ユーザが本システムへの操作を行うために用いる。

出力装置１０は、ディスプレイ、プリンタなどであり、本システムにより生成された林相区画や林相区分図（林分図）を画面表示、印刷等によりユーザに示す等に用いられる。また、林相区画や林相区分に関するデータを他のシステムで利用できるよう、データとして出力してもよい。

林相区画生成部２０は、オルソ画像データ４０に対して階層的に領域分割処理を行って林相区画を生成する。林相区画生成部２０は、低分解能画像生成部７０、初期区画部７２及び後続区画部７４を有する。

低分解能画像生成部７０は、高分解能のオルソ画像データ４０をディグラデーション処理して分解能を低下させた低分解能画像（劣化画像）を生成する画像劣化手段であり、低分解能画像の画素はオルソ画像データ４０における隣接する複数画素からなる領域に対応する。ここでは区別を容易にするために、オルソ画像データ４０の画素を原画素、低分解能画像の画素を拡大画素と表現する。例えば、低分解能画像はオルソ画像データ４０を一様なサンプリング密度でダウンサンプリングして生成することができ、拡大画素の画素値はオルソ画像データ４０にて対応する領域内の複数の原画素を代表する値に設定することができる。例えば、低分解能画像生成部７０は、原画像に対し所定のフィルタを用いた畳み込みにより平滑化処理を行い、拡大画素の画素値を当該拡大画素内の所定位置の原画素での平滑化処理後の画素値に設定することができる。一般的な画像では、注目画素に距離が近い画素の輝度値は注目画素の輝度値と近くなる場合が多く、注目画素から遠ざかるほど、注目画素の輝度値との差が大きくなる場合が多い。このような性質を考慮し、注目画素に近いほど平均値を計算するときの重みを大きくし、遠くなるほど重みを小さくするように２次元ガウス分布を用いて画像の平滑化を行うガウシアンフィルタが知られている。拡大画素の画素値を求める際の平滑化処理に、このガウシアンフィルタを採用することができる。また、より簡単な手法として例えば、拡大画素の画素値はそれに対応する複数の原画素の画素値の平均値、中央値などで定義することもできる。

初期区画部７２は階層的な領域分割処理の最初の階層の処理を行う。具体的には画像のスペクトル特徴量及び幾何特徴量のいずれか、又は両方を区画用特徴量とし、初期区画部７２は低分解能画像生成部７０により生成された低分解能画像を区画用特徴量の類似性に基づいて複数領域に分割し１次の林相区画を生成する。例えば、初期区画部７２は、互いに隣接する複数の拡大画素からなる領域を１次の林相区画とするか否かを、それら拡大画素の画素値のスペクトル特徴量の類似性（スペクトル条件）と、それら拡大画素を結合して得られる領域についての幾何特徴量の類似性（幾何条件）とに基づいて決定する。

後続区画部７４は階層的な領域分割処理の第２の階層以降の処理を行う。具体的には後続区画部７４は、隣接する林相区画を結合して新たな林相区画を生成する逐次区画処理を１回又は複数回行う。逐次区画処理は、互いに隣接する複数の低次の林相区画について、それらのスペクトル特徴量の類似性と、それらを結合して得られる高次の林相区画についての幾何特徴量の類似性とに基づいて、それら低次の林相区画を結合するか否かを決める。

本実施形態では後続区画部７４は逐次区画処理として、初期区画部７２にて生成された低分解能画像の領域分割処理結果である林相区画内における低分解能画像と高分解能画像（原画像）との両方、あるいは高分解能画像のみを用いて領域分割処理を行う。低分解能画像からは、林相区画の結合可否の判断に用いる画素値のスペクトル特徴量として、林相区画内に包含される拡大画素群についてのスペクトル特徴量を求める。一方、高分解能画像からは、林相区画の結合可否の判断に用いる画素値のスペクトル特徴量として、林相区画内に包含される原画素群についてのスペクトル特徴量を求める。なお、後続区画部７４は少なくとも１回の逐次区画処理が、高分解能画像の画素値のスペクトル特徴量、又は高分解能画像及び低分解能画像それぞれの画素値のスペクトル特徴量に基づいて行う領域分割処理である構成とすることもできる。すなわち、後続区画部７４が複数回の逐次区画処理を行う場合に、そのうち一部の回数は低分解能画像の画素値だけからスペクトル特徴量を抽出し、それに基づいて領域分割処理を行ってもよい。

画像領域分割には、例えば領域併合（region merging）に基づく手法を利用することができる。この手法において画像領域分割処理の対象とされる画像に存在する隣接する２つのオブジェクトを統合するかどうかは、統合後に生成される新しいオブジェクトの異質性（heterogeneity）と統合前のオブジェクトの異質性との間の変化を評価することによって決定される。ちなみに、初期区画部７２での処理ではオブジェクトは低分解能画像の画素であり、後続区画部７４での処理ではオブジェクトは既に生成されている林相区画である。

具体的には、領域併合によるオブジェクトの異質性の変化Δｈは、併合前後におけるスペクトル異質性の変化Δｈ_ｐと形状異質性の変化Δｈ_ｔとから次式によって算出される。

ここで、ｗ_ｐはスペクトル異質性の重み、ｗ_ｔは形状異質性の重みである。

併合前後のスペクトル異質性の変化Δｈ_ｐは、対象画像の各バンドにおける併合前後のオブジェクト内の画素値の標準偏差を用いて、次式によって計算される。

ここで、Ｎは画像のバンド数、ｗ_ｉはバンドｉの重み、ｎ_ａｂは併合後の新しいオブジェクトの画素数、ｎ_ａ，ｎ_ｂは併合前の２つのオブジェクトの画素数、σ_ｉ，ａｂは併合後のオブジェクトのバンドｉにおける標準偏差、σ_ｉ，ａ，σ_ｉ，ｂは併合前の２つのオブジェクトのバンドｉにおける標準偏差である。

また、併合前後の形状異質性の変化Δｈ_ｔは、コンパクトネスとスムースネスという２つの基準で次式により定義される。

ここで、Δｈ_ｃは併合前後のコンパクトネスの変化、Δｈ_ｓは併合前後のスムースネスの変化、ｗ_ｃはコンパクトネスの重み、ｗ_ｓはスムースネスの重みである。

オブジェクトのコンパクトネス基準はオブジェクトの周囲長と面積から計算され、一方、スムースネス基準はオブジェクトの周囲長と境界ボックスの直径（長軸）から計算される。具体的にはΔｈ_ｃ，Δｈ_ｓは次式で定義される。

ここで、ｌ_ａｂは併合後のオブジェクトの周囲長、ｌ_ａ，ｌ_ｂは併合前の２つのオブジェクトの周囲長、ｓ_ａｂは併合後のオブジェクトの面積、ｓ_ａ，ｓ_ｂは併合前の２つのオブジェクトの面積、ｂ_ａｂは併合後のオブジェクトの境界ボックスの直径、ｂ_ａ，ｂ_ｂは併合前の２つのオブジェクトの境界ボックスの直径、ｎ_ａｂは併合後の新しいオブジェクトの画素数、ｎ_ａ，ｎ_ｂは併合前の２つのオブジェクトの画素数である。

併合前後のオブジェクトの異質性の変化Δｈが設定されたしきい値を超えない場合、領域の併合処理が実施され、しきい値を上回る場合、領域の併合処理を停止する。設定されたしきい値はスケールパラメータ（scale parameter）と呼ばれ、画像の分割処理によって生成されるオブジェクトの大きさを表している。スケールパラメータが大きいほど、より多くのオブジェクトが併合され、領域分割によって最終的に生成されるオブジェクトのサイズが大きくなる。

（１）式に示すように、初期区画部７２における拡大画素の結合の判断、又は後続区画部７４における低次の林相区画の結合の判断に対するスペクトル特徴量及び幾何特徴量それぞれの寄与比率は重みｗ_ｐ，ｗ_ｃにより調節することができる。ここで、初期区画部７２及び後続区画部７４の一方又は両方は、林相区画の生成判断において、幾何条件を用いずスペクトル特徴量の類似性だけに基づいて行う構成にすることもできる。

テクスチャ解析部２２は、オルソ画像データ４０からテクスチャ情報を抽出するテクスチャ情報抽出手段である。テクスチャ解析部２２は、日向／日陰分布パターン情報抽出部８０、日向／日陰境界パターン情報抽出部８２及びＬＢＰ画像情報抽出部８４を有する。

日向／日陰分布パターン情報抽出部８０はオルソ画像データ４０を二値化処理して日向と日陰とが区別された日向／日陰分布パターンをテクスチャ情報として抽出する。図２は日向／日陰分布パターンの例を示す説明図である。図２（ａ）〜（ｄ）はそれぞれスギ林、ヒノキ林、広葉樹林、非森林の例であり、左右に並ぶ２つの画像のうち左側がオルソ画像データ４０、右側が日向／日陰分布パターンである。当該日向／日陰分布パターンにおいて白領域が日向、黒領域が日陰である。ちなみに、ここでの二値化のしきい値は大津の手法により決定している。

日向／日陰境界パターン情報抽出部８２はオルソ画像データ４０における輝度勾配を二値化処理して日向と日陰との境界領域が現れた日向／日陰境界パターンをテクスチャ情報として抽出する。図３は日向／日陰境界パターンの例を示す説明図である。図３（ａ）〜（ｄ）は図２と同様、それぞれスギ林、ヒノキ林、広葉樹林、非森林の例であり、左右に並ぶ２つの画像のうち左側がオルソ画像データ４０、右側が日向／日陰境界パターンである。当該日向／日陰境界パターンにおいて白領域が輝度勾配がしきい値以上の領域であり、黒領域がしきい値未満の領域である。ここでも、しきい値は大津の手法により決定している。日向／日陰境界パターンからは日向と日陰とが切り替わる空間的な頻度の多寡が読み取れる。

ＬＢＰ画像情報抽出部８４はオルソ画像データ４０から局所二値パターン演算子を用いて得られるＬＢＰ（Local Binary Pattern）画像をテクスチャ情報として抽出する。図４はＬＢＰ画像の例を示す説明図である。図４（ａ）〜（ｄ）は図２及び図３と同様、それぞれスギ林、ヒノキ林、広葉樹林、非森林の例であり、左右に並ぶ２つの画像のうち左側がオルソ画像データ４０、右側がＬＢＰ画像である。ＬＢＰ画像は原画像の詳細な模様構造パターンを反映しており、しかも画像のコントラストの影響を受けにくいという特性を有する。

テクスチャ解析部２２は抽出された日向／日陰分布パターン、日向／日陰境界パターン及びＬＢＰ画像の特徴量として、それぞれの平均値や標準偏差などの統計量を算出する。

スペクトル解析部２４はオルソ画像データ４０からスペクトル情報を抽出するスペクトル情報抽出手段である。本実施形態ではスペクトル解析部２４はＲ，Ｇ，Ｂ，ＮＩＲの４バンドのマルチスペクトル画像に対して次式で示される正規化処理を行い、正規化後の成分Ｒ’，Ｇ’，Ｂ’からなる画像からスペクトル特徴量を抽出する。またスペクトル解析部２４は日向／日陰分布パターン情報に基づいて各林相区画における日向領域のみからスペクトル特徴量を抽出する。

図５及び図６はオルソ画像データ４０に対する正規化処理の例を示す画像であり、図５は原画像、すなわち正規化処理前のオルソ画像、図６は図５に対応するオルソ画像データ４０を正規化処理して得られた画像である。図５及び図６はグレースケールで表示しているが本来、図５はＲＧＢ合成のカラー画像であり、図６はＲ’Ｇ’Ｂ’合成のカラー画像である。例えば、樹種等の違いに応じた近赤外成分の大きさの相違により、原画像と正規化処理画像との間には色合いなどに違いが生じる。図５と図６とにおける原画像と正規化処理画像との相違は、画像内での濃淡の変化が両画像間にて相違していることから読み取れる。

なお、スペクトル解析部２４は正規化処理を行わずにスペクトル特徴量を抽出する構成や、日陰を含んだ各林相区画の全体におけるスペクトル特徴量を抽出する構成とすることもできる。

統計分布作成部２６は林相区分図データ４２として与えられている時期Ｔ１の既知の林相区分図での各林相の領域に対応して、オルソ画像データ４０で与えられる時期Ｔ２の原画像に注目領域を設定し、各林相の注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める。

特徴情報は原画像からテクスチャ解析部２２、スペクトル解析部２４により抽出されるテクスチャ情報、スペクトル情報である。本実施形態では統計分布を求める画像の特徴情報としてテクスチャ情報を用いることに対応して、例えば、原画像について生成した林相区画を要素領域とする。なお、統計分布を求める画像の特徴情報としてスペクトル情報だけを用いるような場合には、画素を要素領域とすることもできる。

例えば、時期Ｔ１の林相区分図におけるスギ林の林相区画の集合に対応して、原画像にスギ林に関する注目領域が設定される。当該注目領域は、時期Ｔ１の林相区分図におけるスギ林の林相区画の集合からなる領域と同一又は近似的に同一の領域である。具体的には、原画像について生成した林相区画のうち、時期Ｔ１のスギ林領域に対応する位置に存在するものの集合をスギ林に関する注目領域として設定することができる。ここで、原画像の或る林相区画が時期Ｔ１の複数の林相の領域に跨がる場合には、例えば、当該林相区画は、時期Ｔ１の林相別領域のうち当該林相区画との重なりが最も大きいものに対応した林相の注目領域に含めることができる。また、重なりの面積の割合に応じて重み付けをして統計分布に反映させることもできる。

統計分布として、画像の特徴情報であるテクスチャ情報やスペクトル情報を表す特徴量を変数とした度数分布を求める。変数は１種類だけであってもよいし複数であってもよい。度数は単純に要素領域の数としてもよいし、林相区画を要素領域とする場合には、林相区画の面積に応じて重み付けをした値を積算してもよい。

変化領域検出部２８は各注目領域での統計分布における外れ値を検出し、当該外れ値を与える要素領域を当該注目領域において時期Ｔ１との間で林相が変化した変化領域として検出する。外れ値は統計において他の値から大きく外れた値である。ここで、既に述べたように通常、林相区分図は数年程度ごとに更新される。つまり、時期Ｔ１と時期Ｔ２との間隔は基本的に数年程度である。林相解析システム２は航空機などによる撮影されたオルソ画像データ４０を用いて広範囲の林相解析を効率的に行うものであり、その広範囲の対象地域にて数年程度の短期間に大規模な林相変化が起こることはあまりない。つまり、通常は、時期Ｔ１と時期Ｔ２との間にて林相が変化する領域は広範囲な対象地域のごく一部であることが期待でき、時期Ｔ２の原画像について設定された各注目領域において、既存の林相区分図が得られている時期Ｔ１との間で林相が変わらない無変化領域が大半を占めることを前提として統計的解析を行うことが可能である。そこで、変化領域検出部２８は或る林相に対応した注目領域について得られた統計分布において、無変化領域について所定の形の確率密度関数を仮定し、当該確率密度関数から大きく外れている標本である要素領域を当該林相とは異なる領域、つまり変化領域であると判定する。

或る林相に対応した注目領域における統計分布は基本的に、対象地域にて当該林相を有する全ての林相区画における特徴情報を反映し、林相が共通する林相区画における特徴情報のばらつきの影響を受けにくい。よって、当該統計分布に基づいて変化領域を判定することで、信頼性が高い変化領域が得られる。

林相区分図生成部３０は、時期Ｔ１の林相区分図に変化領域を反映させて時期Ｔ２の林相区分図を生成する。具体的には、林相区分図生成部３０は、時期Ｔ１の林相区分図における各林相に対応して設定した時期Ｔ２の原画像における注目領域にて、変化領域と判定された林相区画は伐採や風倒被害等による非森林とし、無変化領域である残りの林相区画は当該注目領域に対応する林相に分類し、時期Ｔ２の林相区分図を生成する。ちなみに、時期Ｔ１，Ｔ２の間隔は数年程度の期間である。当該期間は森林の生長からすれば比較的短く、当該期間に例えば、スギ林がヒノキ林に変化するといった林相の変化は基本的にないので、ここでは変化領域は非森林に分類している。なお、各林相についてテクスチャやスペクトルについての特徴情報を例えば、予め記憶装置６に格納しておき、これを用いて変化領域がどの林相の特徴を有しているかを判定し、当該変化領域の林相を決定してもよい。

図７は林相解析システム２における概略のデータフロー図であり、以下、図７に沿って林相解析システム２の処理を説明する。

林相区画生成部２０は時期Ｔ２のオルソ画像データ４０から時期Ｔ２における林相区画データ１００を生成する。図８は林相区画生成部２０による林相区画の生成処理を説明する概略のフロー図である。林相区画生成部２０は低分解能画像生成部７０によりオルソ画像データ４０から低分解能画像を生成する（Ｓ１０）。さらに林相区画生成部２０は低分解能画像を領域分割した林相区画を階層的に領域分割して、内部が一様な林相からなる領域に対応した林相区画を生成する（Ｓ２０）。この領域分割処理Ｓ２０は初期区画部７２及び後続区画部７４により行われる。領域分割処理は低分解能画像から開始され、低分解能画像を構成する拡大画素の画素値を用いた区画処理が最初に少なくとも１回行われ、その後、オルソ画像データ４０を構成する原画素及び拡大画素あるいは原画素のみの画素値を用いた区画処理が少なくとも１回行われる。具体的には、初期区画部７２は低分解能画像を領域分割してスケールが小さい林相区画を生成する（Ｓ２２）。後続区画部７４は低分解能画像と高分解能画像とを併用し、林相区画同士を結合することによりスケールが大きくなった林相区画を生成する（Ｓ２４）。

領域分割処理はスペクトル条件４４及び幾何条件４６を満たすように行われる。その際、スペクトル条件４４と幾何条件４６との比重、つまりそれぞれを領域分割に寄与させる度合は調節することができる。

ここで、領域分割処理の階層ごとに、林相区画を生成する際の条件は異なり得る。すなわち、初回の拡大画素を結合して初期の林相区画を生成する際や、低次の林相区画同士を結合して高次の林相区画を生成する際に、スケールパラメータであるΔｈは生成される林相区画のスケールが徐々に大きくなることを可能とするように設定される。Δｈの増加により、スペクトル条件４４に関しては、結合対象となる複数の拡大画素又は林相区画についてのスペクトル特徴量の類似性の判断基準（Δｈ_ｐ）が領域分割処理が高階になるにつれ緩和される。また、幾何条件４６に関しても幾何特徴量の類似性の判断基準（Δｈ_ｔ）が緩和され、領域分割処理が高階になるにつれより複雑な形状あるいは大きな面積を有する林相区画が許容されるようになり結合が促進される。一方、スペクトル条件４４と幾何条件４６とで総合的に結合が緩和されればよいので、スペクトル条件４４により結合を緩和しつつ、幾何条件４６は結合を或る程度抑制する方向に変化させより単純な形状の林相区画の生成を促す条件に設定することもできる。例えば、そのような調整は（１）式の重みｗ_ｐ，ｗ_ｔを変えることで可能であり、具体的にはｗ_ｐを低下させ、一方、ｗ_ｔを増加させて、Δｈ_ｐがΔｈに寄与しにくくし、Δｈ_ｔがΔｈに寄与しやすくすることで実現できる。

さて、既に述べたように、高分解能画像では森林の日向と日陰の違いにより同じ林相内でも画素値のばらつきが大きいので、従来、隣接している林相を区別して正確に区画するのが困難であるという問題があった。これに対し、林相解析システム２は原画像を低分解能化し、階層的な領域分割処理のうち少なくとも最初の１回を低分解能画像にて行う。画像を低分解能化すると同じ林相内での画素値のばらつきが平滑化・緩和されるので、低分解能画像に基づく区画処理は日向と日陰とでの画素値の違いの影響を受けにくくなり、隣接する異なる林相を区別し易い。なお、低分解能での区画処理を階層的な複数回の区画処理のうちの早い段階で行う方が、当該効果は高くなることが期待できる。

一方、階層的な区画処理の中で、低分解能画像と高分解能画像との両方、あるいは高分解能画像のみを用いて区画処理を行うことで、最終的に得られる林相区画に高分解能画像の情報を反映させることができる。

図９〜図１１は林相区画生成部２０の処理の例を示す画像であり、それぞれ同じ森林の空中写真画像に林相区画の境界を黒線で描いている。図９は初期区画部７２による処理結果であり、低分解能画像を利用して生成された初期の林相区画を示している。図１０及び図１１は後続区画部７４による処理結果であり、図１１の画像は図１０の画像より高階の領域分割処理の結果を示している。すなわち、図９、図１０、図１１の順に階層的に領域分割が行われる。ここで、上述したように林相区画の生成は低分解能画像と高分解能画像とを併用して行われる。具体的には、図９〜図１１に示す３階層が林相区画を生成する領域分割処理の全階層である場合には、図９の林相区画は低分解能画像を利用して生成され、図１０及び図１１の林相区画は、先行する領域分割の処理結果である林相区画内における低分解能画像と高分解能画像との両方、あるいは高分解能画像のみを利用して生成される。そして、図１１の林相区画が最終的に抽出された林相区画となる。

上述した林相区画生成部２０の処理により対象地域から最終的に抽出された林相区画のデータ１００はテクスチャ解析部２２及びスペクトル解析部２４、並びに林相区分図生成部３０にて利用される。

さて、説明を図７に戻す。テクスチャ解析部２２及びスペクトル解析部２４はオルソ画像データ４０及び林相区画データ１００とから、時期Ｔ２における林相区画ごとに画像の特徴情報データ１０２を抽出する。図１２はテクスチャ解析部２２及びスペクトル解析部２４による特徴情報データ１０２の抽出処理を説明する概略のフロー図である。

林相解析システム２はテクスチャ解析部２２によりオルソ画像データ４０からテクスチャ情報を抽出する（Ｓ４０，Ｓ５０，Ｓ６０）。具体的には、テクスチャ解析部２２は日向／日陰分布パターン情報抽出部８０により日向／日陰分布パターンを抽出し（Ｓ４０）、日向／日陰境界パターン情報抽出部８２により日向／日陰境界パターンを抽出し（Ｓ５０）、またＬＢＰ画像情報抽出部８４により、原画像における詳細な模様構造パターンが現れるＬＢＰ画像を抽出し（Ｓ６０）、テクスチャ情報を求める。

さらに林相解析システム２はスペクトル解析部２４によりオルソ画像データ４０を正規化処理し、正規化後の画像からスペクトル特徴情報を抽出する（Ｓ７０）。さらにスペクトル解析部２４は日向／日陰分布パターン情報抽出部８０により抽出した日向／日陰分布パターンを用いて、森林のスペクトルの特徴が良好に現れ得る日向部分のみからスペクトル特徴情報を抽出する（Ｓ８０）。

そしてテクスチャ解析部２２及びスペクトル解析部２４は林相区画生成部２０により抽出された林相区画ごとに画像の特徴情報を取得し（Ｓ９０）、特徴情報データ１０２を生成する。この特徴情報データ１０２は統計分布作成部２６にて利用される。

再び図７に戻る。統計分布作成部２６は時期Ｔ１の林相区分図データ４２と時期Ｔ２の林相区画ごとの画像の特徴情報データ１０２とから、時期Ｔ２の原画像にて林相ごとに設定される注目領域それぞれについて画像の特徴情報の統計分布を求める。変化領域検出部２８は統計分布作成部２６から出力される統計分布データ１０４に基づいて、各注目領域に含まれる変化領域を示す変化領域データ１０６を生成する。

図１３は１つの林相に対応した注目領域に関する統計分布の模式図であり、横軸が特徴情報に関する特徴量Ｘで表される変数に対応し、縦軸が度数Ｙに対応する。ここでは図示の都合上、統計分布１０８は変数が１つである１次元の分布としているが、統計分布作成部２６は変数が複数である多次元の統計分布を求めてもよい。例えば、変数はユーザにより、注目領域に対応する林相と、当該注目領域にて時期Ｔ１，Ｔ２間に生じることが想定される森林変化、具体的には伐採や風倒被害等に対応する林相とを弁別するのに好適な単一の特徴量又は複数の特徴量の組み合わせに設定される。上述したテクスチャ解析部２２及びスペクトル解析部２４により算出される各特徴情報は当該弁別に好適であり、統計分布の変数はそれら複数種類の特徴情報のうちの一部又は全部で定義することができる。

林相が共通する要素領域（林相区画、画素等）についての特徴量の統計分布は、基本的には平均値の付近に集積してピークを形成するような分布となり、例えば、正規分布を仮定することができる。既に述べたように、対象地域において時期Ｔ１から林相が変化する変化領域は、変化しない領域に比べて少数であるとすることができる。よって、各注目領域についての統計分布１０８は当該注目領域に対応する林相の特徴量の分布に近似したものとなる。つまり当該林相の要素領域は統計分布１０８において最大のピークの近傍の特徴量を有する確率が高い。一方、統計分布の変数とする特徴量Ｘは上述のように林相の変化の有無を弁別容易とするものに設定され、林相が変化した要素領域は統計分布１０８においてピークから離れた範囲（例えば図１３における範囲Ｒ_Ｅ）に属する特徴量を有することが期待できる。そこで、変化領域検出部２８は或る林相に対応した注目領域について得られた統計分布１０８において外れ値となる要素領域を検出し、当該要素領域を示す情報を変化領域データ１０６として出力する。例えば、林相区画生成部２０は生成した林相区画を識別する識別子（ＩＤ）を定義し、林相区画データ１００及び特徴情報データ１０２はそれぞれ当該ＩＤに対応付けられた林相区画の情報、特徴情報とされる。変化領域検出部２８は変化領域データ１０６として例えば、当該ＩＤを抽出することができる。

変化領域検出部２８は例えば、或る林相に対応する注目領域の統計分布が当該林相の特徴情報について想定される分布に従っているかを検証する仮説検定を用い、仮説を棄却する原因となっている標本を探索して、それを外れ値とする。

例えば、各林相の特徴量の統計分布が正規分布である場合に、カイ２乗（χ^２）検定を用いることができる。変化領域検出部２８は、林相の特徴情報が従うべき正規分布の平均値の推定値として注目領域の統計情報の平均値を用いてχ^２値を計算する。χ^２値がχ^２分布において所定の有意水準（例えば、１％）での棄却域にある場合は、注目領域にて観測された特徴情報からなる標本の集合うち分散が最大のものを除去した残りの集合についてχ^２値を計算し直す。この作業をχ^２値が採択域に入るまで反復する。変化領域検出部２８は、この処理にて元の標本の集合から除去された標本を外れ値とし、変化領域として検出する。

なお、外れ値の検出は例えば、スミルノフ・グラブス検定などの他の周知の検定手法を用いて行うこともできる。

林相区分図生成部３０は林相区画データ１００で与えられる各林相区画について、時期Ｔ１の林相区分図データ４２に基づいていずれの林相の注目領域に属するかを判定し、さらに変化領域データ１０６に基づいて当該注目領域における変化領域であるか否かを判定する。これにより時期Ｔ２の各林相区画の林相が決定され、時期Ｔ２についての林相区分図（林相区分図データ１１０）が生成される。

なお、テクスチャ情報は上述の３種類以外のものを用いることもできる。例えば、ウェーブレット解析により抽出した画像テクスチャを用いることができる。具体的には、ウェーブレット解析で抽出した高周波成分の情報をテクスチャ情報として画像のスペクトル特徴情報と併用して林相を判定することができる。

本実施形態ではオルソ画像データ４０は航空写真としたが、本発明は、例えば衛星写真などの他の高分解能画像を用いた林相解析にも適用可能である。また、本実施形態ではオルソ画像を例に説明したが、必ずしもオルソ画像を用いる必要は無く、上空から撮影した画像を適宜用いることができる。

［第２の実施形態］
第２の実施形態の林相解析システムにおいて上記第１の実施形態と同じ構成要素は同一の符号を付して説明を省略する。以下、第２の実施形態の林相解析システム２００について第１の実施形態の林相解析システム２との相違点を説明する。

図１４は第２の実施形態である林相解析システム２００の概略の構成を示すブロック図である。林相解析システム２００と林相解析システム２との相違点は林相区画をどのような条件に基づいて行うかという点にある。第１の実施形態では林相区画生成部２０はスペクトル特徴量についてのスペクトル条件４４及び幾何特徴量についての幾何条件４６を満たすように領域分割処理を行って林相区画を生成している。これに対し、第２の実施形態の林相解析システム２００では、領域併合前後のオブジェクト、つまり、劣化画像の画素又は林相区画の色、形状、テクスチャ及び面積を区画用特徴量とし、これらについての色条件２０２、形状条件２０４、テクスチャ条件２０６及び面積条件２０８を満たすように領域分割処理を行って林相区画を生成する。

色条件２０２、形状条件２０４、テクスチャ条件２０６及び面積条件２０８は、オルソ画像（オルソ画像データ４０）を領域分割して生成されるオブジェクト領域である林相区画が満たすべき条件である。色条件２０２は林相区画におけるオルソ画像の画素値に関する条件を規定する。ちなみに「色条件」における「色」はオルソ画像の画素値を意味し、ＲＧＢの情報からなるスペクトル情報よりも広義で用いている。具体的には、ここでの「色」は第１の実施形態のスペクトル特徴量に相当し、例えば、ＮＩＲの情報を含むマルチスペクトル画像の画素値もここでの色の概念に含まれ得る。形状条件２０４は林相区画の形状に関する条件を規定する。テクスチャ条件２０６は林相区画の画像テクスチャに関する条件を規定する。画像テクスチャはテクスチャ特徴量で評価される。面積条件２０８は林相区画の面積に関する条件を規定する。なお、林相区画の形状、面積は林相区画の幾何特徴量である。

初期区画部７２は、互いに隣接する複数の画素からなる領域を１次の林相区画とするか否かを、それら画素の画素値の類似性、領域の画像テクスチャの類似性、及び領域の形状や面積などの幾何特徴量の類似性に基づいて決定する。ここで、結合される領域間や結合前後の領域間での画素値の類似性は、画素値の平均値などの統計量などの特徴量に基づいて判断することができる。

後続区画部７４は第１の実施形態と同様、逐次区画処理を行う。本実施形態では逐次区画処理は、互いに隣接する複数の低次の林相区画について、それらにおける画素値の類似性に関する色条件２０２と、それらを結合して得られる高次の林相区画についての形状条件２０４と、テクスチャ条件２０６と、面積条件２０８とに基づいて、それら低次の林相区画を結合するか否かを決める。

画像領域分割には第１の実施形態と同様、領域併合に基づく手法を利用する。本実施形態では、領域併合によるオブジェクトの異質性の変化Δｈは、併合前後における色の異質性の変化Δｈ_colorと、形状の異質性の変化Δｈ_shapeと、テクスチャの異質性の変化Δｈ_textureと、面積の異質性の変化Δｈ_aeraとから次式によって算出される。

ここで、ｗ_colorは色の異質性の重み、ｗ_shapeは形状の異質性の重み、ｗ_textureはテクスチャの異質性の重み、ｗ_areaは面積の異質性の重みである。

併合前後の色の異質性の変化Δｈ_color、形状の異質性の変化Δｈ_shape、テクスチャの異質性の変化Δｈ_texture、面積の異質性の変化Δｈ_aeraは、例えば、それぞれ次式によって計算される。

ここで、ｎ_ａｂは併合後の新しいオブジェクトの画素数、ｎ_ａ，ｎ_ｂは併合前の２つのオブジェクトの画素数である。また、ＣＩ_ａｂは併合後のオブジェクトの色情報指数（例えば、画素値）、ＣＩ_ａ，ＣＩ_ｂは併合前の２つのオブジェクトの色情報指数、ＳＩ_ａｂは併合後のオブジェクトの形状情報指数（例えば、スムースネス）、ＳＩ_ａ，ＳＩ_ｂは併合前の２つのオブジェクトの形状情報指数、ＴＩ_ａｂは併合後のオブジェクトのテクスチャ情報指数（例えば、画素値の標準偏差値）、ＴＩ_ａ，ＴＩ_ｂは併合前の２つのオブジェクトのテクスチャ情報指数、ＡＩ_ａｂは併合後のオブジェクトの面積情報指数（例えば、面積値）、ＡＩ_ａ，ＡＩ_ｂは併合前の２つのオブジェクトの面積情報指数である。

（６）式に示すように、初期区画部７２における画素の結合の判断、又は後続区画部７４における低次の林相区画の結合の判断に対する色条件、形状条件、テクスチャ条件及び面積条件それぞれの寄与比率は重みｗ_color，ｗ_shape，ｗ_texture，ｗ_areaにより調節することができる。ここで、初期区画部７２及び後続区画部７４の一方又は両方は、林相区画の生成判断において、色条件、形状条件、テクスチャ条件及び面積条件の何れか１つもしくは複数の組み合わせに基づいて行う構成にすることもできる。

２，２００ 林相解析システム、４ 演算処理装置、６ 記憶装置、８ 入力装置、１０ 出力装置、２０ 林相区画生成部、２２ テクスチャ解析部、２４ スペクトル解析部、２６ 統計分布作成部、２８ 変化領域検出部、３０ 林相区分図生成部、４０ オルソ画像データ、４２ 時期Ｔ１の林相区分図データ、４４ スペクトル条件、４６ 幾何条件、７０ 低分解能画像生成部、７２ 初期区画部、７４ 後続区画部、８０ 日向／日陰分布パターン情報抽出部、８２ 日向／日陰境界パターン情報抽出部、８４ ＬＢＰ画像情報抽出部、１００ 林相区画データ、１０２ 特徴情報データ、１０４ 統計分布データ、１０６ 変化領域データ、１１０ 時期Ｔ２の林相区分図データ、２０２ 色条件、２０４ 形状条件、２０６ テクスチャ条件、２０８ 面積条件。

Claims (9)

1. 森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成手段と、
   前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出手段と、
   を有することを特徴とする林相解析装置。
2. 請求項１に記載の林相解析装置において、
   前記第１の時期の前記林相区分図に前記変化領域を反映させて前記第２の時期の林相区分図を生成する林相区分図生成手段を有すること、を特徴とする林相解析装置。
3. 請求項１又は請求項２に記載の林相解析装置において、
   前記原画像の１又は複数種類の区画用特徴量の類似性に基づいて前記対象地域を複数領域に分割して前記森林の林相区画を生成し出力する林相区画生成手段を有し、
   前記統計分布作成手段及び前記変化領域検出手段における前記要素領域は前記林相区画であること、
   を特徴とする林相解析装置。
4. 請求項３に記載の林相解析装置において、
   前記原画像から、当該原画像の隣接する複数画素からなる領域を画素として分解能を低下させた劣化画像を生成する画像劣化手段を有し、
   前記林相区画生成手段は、
   前記劣化画像を、その１又は複数種類の前記区画用特徴量の類似性に基づいて複数領域に分割し林相区画を生成する初期区画手段と、
   隣接する林相区画をそれらの１又は複数種類の前記区画用特徴量の類似性に基づいて結合して新たな林相区画を生成する逐次区画処理を１回又は複数回行い、前記対象地域における前記森林の林相区画を出力する後続区画手段と、
   を有し、
   前記後続区画手段は少なくとも１回の前記逐次区画処理を前記原画像の前記区画用特徴量、又は前記原画像及び前記劣化画像それぞれの前記区画用特徴量に基づいて行うこと、
   を特徴とする林相解析装置。
5. 請求項３又は請求項４に記載の林相解析装置において、
   前記原画像からテクスチャ情報を抽出するテクスチャ情報抽出手段を有し、
   前記統計分布作成手段は、前記林相区画生成手段から出力される前記各林相区画から抽出された前記テクスチャ情報を前記特徴情報として前記統計分布を求めること、
   を特徴とする林相解析装置。
6. 請求項５に記載の林相解析装置において、
   前記原画像からスペクトル情報を抽出するスペクトル情報抽出手段を有し、
   前記統計分布作成手段は、前記林相区画生成手段から出力される前記各林相区画から抽出された前記テクスチャ情報及び前記スペクトル情報を前記特徴情報として前記統計分布を求めること、
   を特徴とする林相解析装置。
7. 請求項３又は請求項４に記載の林相解析装置において、
   前記原画像からスペクトル情報を抽出するスペクトル情報抽出手段を有し、
   前記統計分布作成手段は、前記林相区画生成手段から出力される前記各林相区画から抽出された前記スペクトル情報を前記特徴情報として前記統計分布を求めること、
   を特徴とする林相解析装置。
8. 森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成ステップと、
   前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出ステップと、
   を有することを特徴とする林相解析方法。
9. コンピュータに、森林を含む対象地域における林相解析を行わせるためのプログラムであって、当該コンピュータを、
   森林を含む対象地域の第１の時期における既知の林相区分図での各林相の領域に対応して、前記対象地域を第２の時期に上空から撮影した原画像に注目領域を設定し、前記各林相の前記注目領域にて所定の要素領域ごとに抽出される画像の特徴情報について統計分布を求める統計分布作成手段、及び、
   前記各注目領域での前記統計分布における外れ値を検出し、当該外れ値を与える前記要素領域を当該注目領域において前記第１の時期との間で前記林相が変化した変化領域として検出する変化領域検出手段、
   として機能させることを特徴とするプログラム。

Images