JP6200826B2

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

Info

Publication number: JP6200826B2
Application number: JP2014024104A
Authority: JP
Inventors: 秀樹島村; 林朱; スバスチャタクリ
Original assignee: Pasco Corp
Current assignee: Pasco Corp
Priority date: 2014-02-12
Filing date: 2014-02-12
Publication date: 2017-09-20
Anticipated expiration: 2034-02-12
Also published as: JP2015152340A

Description

本発明は、航空画像及び航空レーザ計測データから森林の林相区画における林相を判別する林相解析装置、林相解析方法及びプログラムに関する。

林相は樹種・樹齢、樹冠や木の生育状態などによる森林の様子・形態であり、林相区画は林相によって区分された森林領域である。一般に、林相区分図は林相区画における林相のタイプ（林相区分）を判別することで作成される。

ここで、航空機などによる上空からのリモートセンシングは、広範囲での地上の状況把握が可能であり、林相区分図の作成などの林相の解析に利用されている。具体的には、従来、林相の解析には主として航空写真等の上空からの高分解能画像データが利用されている。そして当該データを利用した林相区画の抽出や林相区分図の作成は、人間が実体視鏡やデジタル図化機等を用いて目視判読で行っている。また、人間による広範囲にわたる林相判読は膨大な労力・時間を要するため、航空写真等の画像の色やテクスチャの情報を用いて自動的に林相解析を行う手法の研究が進められている。しかし、航空写真を用いた林相解析は、撮影の時期、撮影の方向・角度、及び太陽の方位角・高度角などにより画像が影響を受け、またそれを避けようとすると撮影機会が制限されるという問題を有する。また、例えば、樹種が異なる林相を色で区別することが難しい場合があることや、同じ樹種であっても樹齢や生育状況の異なる林相を、高さ情報を含まない航空写真から区別することは難しいことなどから、航空写真に基づいて正確な林相区画を生成することは必ずしも容易ではない。

ここで他のリモートセンシング技術として航空レーザ計測がある。航空レーザ計測で取得した３次元点群データは森林の地形や樹高の計測、立木密度や材積の推定などに幅広く利用されている。航空レーザ計測データからは樹木の高さ情報を得ることができるので、航空写真と航空レーザ計測データとを併用して林相解析を行うことも研究されている。しかし従来提案されている林相解析では、航空レーザ計測データに含まれ得る情報のうち高さや反射強度など単純な情報しか利用しないものが主である。

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

航空レーザ計測は地物の３次元構造情報を直接計測によって取得し、またレーザを自ら照射する能動的なセンシングであるので、航空写真に比べて地物の３次元構造情報データ収集時に関する条件が緩やかである。そのため、航空写真と航空レーザ計測とを併用することで、より安定かつ正確に林相区分図の作成などの林相解析を行える可能性がある。しかし、航空写真や衛星画像などの空中撮影画像（アナログ写真、デジタル写真を含む）及び航空レーザ計測データを併用した林相解析において、空中撮影画像及び航空レーザ計測データそれぞれからどのような情報を抽出し、どのように活用すると好適な結果が得られるかは解決すべき課題であった。

本発明は、空中撮影画像及び航空レーザ計測データを併用した林相解析により、好適な林相区分図の自動生成を可能とすることを目的とする。

（１）本発明に係る林相解析装置は、森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成手段と、前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成手段と、前記空中撮影画像及び前記航空レーザ計測データに基づき、予め定められた１又は複数種類の林相特徴量を求める林相特徴量抽出手段と、前記林相区画ごとの前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別手段と、を有する。

（２）上記（１）に記載する林相解析装置において、前記区画用特徴量は、レーザの照射パルスに対する反射パルスの個数に応じて定まる前記各計測地点での反射パルス指標、レーザの反射強度及び樹高のうちの少なくとも１つを含む構成とすることができる。

（３）上記（１）及び（２）に記載する林相解析装置において、前記林相特徴量は、レーザの照射パルスに対する反射パルスの個数に応じて定まる前記各計測地点での反射パルス指標、レーザの反射強度、当該レーザの反射強度を画素値とした前記対象地域での画像を二値化したパターンから抽出されたテクスチャ特徴量、及び樹高のうちの少なくとも１つと、それぞれ前記空中撮影画像から抽出されるテクスチャ情報及びスペクトル情報のうちの少なくとも１つとを含む構成とすることができる。

（４）本発明に係る林相解析方法は、森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成ステップと、前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成ステップと、前記空中撮影画像及び前記航空レーザ計測データに基づき、予め定められた１又は複数種類の林相特徴量を求める林相特徴量抽出ステップと、前記林相区画ごとの前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別ステップと、を有する。

（５）本発明に係るプログラムは、コンピュータに林相解析を行わせるためのプログラムであって、当該コンピュータを、森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成手段、前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成手段、前記空中撮影画像及び前記航空レーザ計測データに基づき、予め定められた１又は複数種類の林相特徴量を求める林相特徴量抽出手段、及び、前記林相区画ごとの前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別手段、として機能させる。

本発明によれば、空中撮影画像及び航空レーザ計測データを併用した林相解析により、好適な林相区分図の自動生成が可能となる。

本発明の実施形態である林相解析システムの概略の構成を示すブロック図である。本発明の実施形態である林相解析システムにおける処理及びデータの概略のフロー図である。航空レーザ計測データから得られる特徴量の例を説明する模式図である。対象地域の航空写真（オルソ画像）の一例である。図４に示す対象地域のＤＣＨＭ（デジタル樹冠高モデル：Digital Canopy Height Model）画像の一例である。図４に示す対象地域のＲＩ（反射強度：Reflection Intensity）画像の一例である。図４に示す対象地域のＢＲＩ（二値化反射強度：Binary Reflection Intensity）画像の一例である。図４に示す対象地域のＦＰＲ（ファーストパルス比：First Pulse Ratio）画像の一例である。図４に示す対象地域のＩＰＲ（中間パルス比：Intermediate Pulse Ratio）画像の一例である。林相区画生成部の概略の処理フロー図である。林相区画生成部により生成された林相区画の一例を示す画像である。画像ＢＩ（二値画像：Binary Image）の例を示す説明図である。画像ＢＧＩ（二値化勾配画像：Binary Gradient Image）の例を示す説明図である。ＬＢＰ（局所二値パターン：Local Binary Pattern）画像の例を示す説明図である。図４に示す対象地域での最近隣法による林相分類で作成された林相区分図である。

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

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

演算処理装置４は、コンピュータのＣＰＵ（Central Processing Unit）からなり、後述するレーザ計測データ解析部２０、空中撮影画像解析部２２、林相区画生成部２４及び林相区分決定部２６として機能する。

記憶装置６はＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等のメモリ装置である。記憶装置６は演算処理装置４にて実行される各種のプログラムや、本システムの処理に必要な各種データなどを記憶し、演算処理装置４との間でこれらの情報を入出力する。例えば、記憶装置６には、航空レーザ計測データ３０、オルソ画像データ３２並びに林相区画に関する色条件３４、形状条件３６、テクスチャ条件３８及び面積条件４０が予め格納され、また演算処理装置４により生成された林相区画データ４２や林相区分図データ４４が格納される。

航空レーザ計測データ３０は、例えば、航空機やヘリコプターなどに搭載されたレーザ計測システムを用いて取得される。レーザ計測システムはレーザスキャナ及びＧＰＳ／ＩＭＵ（Global Positioning System/Inertial Measurement Unit）を含んで構成される。レーザスキャナは上空から地上へ向けてレーザパルスを掃射し、その反射パルスを受信する。航空レーザ計測データ３０の取得に用いるレーザスキャナは近赤外レーザパルスを照射し、また１つの照射パルスに対する所定数（例えば、４点とする装置が多い）の反射パルスを記録できるものを用いる。レーザスキャナによりレーザパルスの照射方向、パルスの発射から受信までの時間差が得られ、一方、ＧＰＳ／ＩＭＵにより航空機の位置・姿勢が得られ、これらのデータからレーザパルスの反射点の座標が算出される。航空レーザ計測データ３０には例えば、レーザパルスの反射強度、レーザパルスごとの各反射点の座標が含まれる。

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

色条件３４、形状条件３６、テクスチャ条件３８及び面積条件４０は、レーザ計測データ解析部２０により生成される特徴量画像を領域分割して生成されるオブジェクト領域である林相区画が満たすべき条件である。色条件３４は林相区画における特徴量画像の画素値に関する条件を規定する。ちなみに「色条件」における「色」は特徴量画像の画素値を意味し、オルソ画像データから得られＲＧＢの情報からなるスペクトル情報よりも広義で用いている。例えば、ＮＩＲの情報を含むマルチスペクトル画像の画素値や、航空レーザ計測データから得られるＤＣＨＭ画像やＲＩ画像の画素値も、ここでの色の概念に含まれ得る。形状条件３６は林相区画の形状に関する条件を規定する。テクスチャ条件３８は林相区画の画像テクスチャに関する条件を規定する。面積条件４０は林相区画の面積に関する条件を規定する。

林相区画データ４２は、森林を含む対象地域における林相区画を表すデータであり、林相区画生成部２４により生成され、林相区分決定部２６における林相区分図の生成処理にて利用される。

林相区分図データ４４は、本システムにより得られる林相区分図のデータであり、林相区分決定部２６により生成される。

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

出力装置１０は、ディスプレイ、プリンタなどであり、本システムにより生成された林相区分図を画面表示、印刷等によりユーザに示す等に用いられる。

レーザ計測データ解析部２０は、記憶装置６に格納される航空レーザ計測データ３０に基づき、予め定められた１又は複数種類の特徴量を、解析対象とする地域（対象地域）の各計測地点について算出し、計測地点間にて内挿処理を行い各画素での特徴量を定義し、対象地域について当該特徴量を画素値とする特徴量画像を生成する。例えば、レーザ計測データ解析部２０には特徴量画像を生成する特徴量の種類ごとに画像生成部が設けられる。本実施形態ではレーザ計測データ解析部２０がＤＣＨＭ画像生成部５０、ＲＩ画像生成部５２、ＢＲＩ画像生成部５４及びＰＲ画像生成部５６を備えた例を示している。これら各画像生成部については後述する。

空中撮影画像解析部２２は、記憶装置６に格納されるオルソ画像データ３２に基づき、予め定められたテクスチャ情報及びスペクトル情報のいずれか又は両方を例えば、オルソ画像データ３２の画素ごとに算出する。本システムでは空中撮影画像解析部２２はテクスチャ情報を求めるテクスチャ解析部６０と、スペクトル情報を求めるスペクトル解析部６２とを備える。これら各解析部については後述する。

林相区画生成部２４は、レーザ計測データ解析部２０にて生成される特徴量画像のうち、林相区画に用いる特徴量（区画用特徴量）に関するものについて、当該特徴量画像上にて特徴量の類似性に基づいて対象地域を複数領域に分割して林相区画を生成する。例えば、林相区画生成部２４は初期区画部７０及び後続区画部７２を有し、特徴量画像に対して階層的に領域分割処理を行って林相区画（林相区画データ４２）を生成する。

林相区分決定部２６は、対象地域に設定された林相区画ごとの特徴量に基づいて、当該区画における林相（林相区分）を判別する林相判別手段である。林相区分決定部２６は、林相区画データ４２から林相区画を取得し、当該林相区画の林相を判別する。

図２は林相解析システム２における処理及びデータの概略のフロー図である。特徴量画像生成処理Ｓ１０は航空レーザ計測データ３０に基づき、区画用特徴量についての特徴量画像を生成する。なお、以下では航空レーザ計測データ３０のみを用いる場合を例に説明するが、オルソ画像データ３２を併用して特徴量画像を生成してもよい。

区画用特徴量には、その値や値の分布範囲に林相に応じて差異が生じる量が用いられる。図３は航空レーザ計測データ３０から抽出される特徴量の例を説明する模式図であり、森林を含む或る対象地域を秋に観測した例を示している。同図には複数の特徴量それぞれについて観測により得られた数値が示されており、例えば、“Ｄ”に対応する数値は広葉樹林での特徴量の値であり、同様に、“Ｈ”，“Ｓ”，“ＮＦ”の数値はヒノキ林、スギ林、非森林領域での値である。

デジタル樹冠高モデル（Digital Canopy Height Model：ＤＣＨＭ）は数値表層モデル（Digital Surface Model：ＤＳＭ）に含まれる数値標高モデル（Digital Terrain Model：ＤＴＭ）の影響を除去し正規化したＤＳＭである。ちなみに、航空レーザ計測データ３０からＤＳＭ及びＤＴＭを生成し、ＤＳＭからＤＴＭを減算してＤＣＨＭが生成される。ＤＣＨＭが森林から想定される或るしきい値（例えば、図３の例では１５ｃｍ程度）より低い場所は、水域、草地、更地などの非森林領域（ＮＦ）とすることができる。また、建物の領域ではＤＣＨＭから得られる高さが建物に想定されるしきい値以上で、高さの分散が森林に比べて小さくなる場合がある。よって、ＤＣＨＭは森林領域と非森林領域である建物領域とを区別し得る。

森林地域ではレーザの反射パルスのうちファーストパルスに基づいてＤＳＭを生成し、ラストパルスに基づいてＤＴＭを生成し、ＤＳＭからＤＴＭを減算してＤＣＨＭが生成される。ＤＣＨＭは樹高を表すデータとして用いられる。図３には示していないが、樹高は森林の生育状況、樹齢といった林相情報を反映する。

レーザパルスの反射強度（Reflection Intensity：ＲＩ）は、樹木の反射断面積や吸収率を反映する。特に、近赤外領域の波長を有する光に対する反射率の相違は、従来より植生のリモートセンシングに利用されている。図３に示す観測結果では、ヒノキ林（Ｈ）及び広葉樹林（Ｄ）がスギ林（Ｓ）より高い値を示している。ちなみに反射強度はファーストパルスの強度に基づいて定義している。なお、計測に用いられるレーザスキャナの照射パルスの強度が一定であれば、反射強度は反射パルスの強度の絶対値で表すことができ、一方、照射パルスの強度が変わり得る場合には、反射強度は照射パルスの強度で正規化した相対値で表すのが好適である。

図３には、反射パルス指標の例として、トータルパルス（Total Pulse：ＴＰ）、ファーストパルス比（First Pulse Ratio：ＦＰＲ）、及び中間パルス比（Intermediate Pulse Ratio：ＩＰＲ）を示している。ここで、レーザの照射パルスに対する反射パルスの個数に応じて定まる値を反射パルス指標と総称している。その指標値は各計測地点にて定義される。マルチパルス機能を有したレーザスキャナにより取得した航空レーザ計測データに基づく反射パルス指標は森林の内部の構造情報を反映し得る。また、レーザ計測システムとして、フルウェーブ（Full-Waveform）計測に対応したものを用いてもよい。

ＴＰは各照射パルスに対する反射パルスの総数に基づく指標であり、例えば、地表の計測地点の単位面積当たりでの反射パルス数とすることができる。なお、上述のようにマルチパルス機能を有したレーザスキャナにより取得したＴＰは基本的に、単位面積の地表への照射パルス数とは必ずしも一致しない。

図３に示すように、ヒノキ林（Ｈ）、スギ林（Ｓ）及び広葉樹林（Ｄ）におけるＴＰは、広葉樹林、スギ林、ヒノキ林の順に小さくなる結果が観測された。

ＦＰＲは、計測地点の例えば単位面積からの反射パルス総数（つまりＴＰ）に対するファーストパルスの数の比の値である。ちなみにファーストパルスは照射パルスに対して最初に検知される反射パルスである。例えば、単位面積に２発のレーザパルスが照射され、１発に対しては２つの反射パルスが受信され、他の１発に対しては３つの反射パルスが受信された場合、ＴＰは５、ファーストパルスの数は２となるので、ＦＰＲは２／５、つまり４０％である。

ＩＰＲは、計測地点の例えば単位面積からの反射パルス総数（つまりＴＰ）に対する中間パルスの数の比の値である。ちなみに中間パルスは照射パルスに対する反射パルスのうち最初に検知されるファーストパルス及び最後に検知されるラストパルスを除いた残りのパルスであり、１つの照射パルスに対する反射パルスの数が１又は２の場合は中間パルスの数は０となる。例えば、単位面積に２発のレーザパルスが照射され、１発に対しては２つの反射パルスが受信され、他の１発に対しては３つの反射パルスが受信された場合、ＴＰは５、中間パルスの数は１となるので、ＩＰＲは１／５、つまり２０％である。

図３の観測結果においてＦＰＲは基本的にＴＰに反比例する傾向を示している。またＩＰＲは基本的にＴＰと正の相関を有するが、樹種間における差異の程度はＴＰとＩＰＲとで異なる。なお、ＦＰＲ、ＩＰＲは無次元量であり、レーザパルスの照射密度の影響を受けにくい。

図２に戻り特徴量画像生成処理Ｓ１０についての説明を続ける。本実施形態では、例えば区画用特徴量として反射パルス指標、反射強度及び樹高を利用する。反射パルス指標に関してはＦＰＲ（又はＩＰＲ）を特徴量として利用する。これに対応して、ＤＣＨＭ画像生成部５０、ＲＩ画像生成部５２及びＰＲ画像生成部５６が特徴量画像生成処理Ｓ１０を行う特徴量画像生成手段として機能する。なお、特徴量画像として、オルソ画像データを解析して求めたテクスチャ情報やスペクトル情報を併用することも可能である。

具体的には、ＤＣＨＭ画像生成部５０は、航空レーザ計測データ３０からＤＣＨＭを生成し、ＤＣＨＭ画像を生成する。ＲＩ画像生成部５２は航空レーザ計測データ３０から反射強度を取得し、ＲＩ画像を生成する。また、ＰＲ（パルス比：Pulse Ratio）画像生成部５６は航空レーザ計測データ３０からパルス比ＰＲとしてＦＰＲ又はＩＰＲを算出し、ＦＰＲ画像又はＩＰＲ画像を生成する。

生成されたＤＣＨＭ画像、ＲＩ画像及び、ＦＰＲ画像又はＩＰＲ画像は特徴量画像８０として林相区画生成処理Ｓ１２にて利用される。

なお、ＲＩ画像を二値化処理して生成するＢＲＩ（Binary Reflection Intensity）画像から得られるテクスチャ特徴量も区画用特徴量として利用することができ、その場合にはＢＲＩ画像生成部５４は特徴量画像生成手段として機能し、ＢＲＩ画像は特徴量画像８０に含められる。ここで、二値化のしきい値は例えば大津の手法により決定することができる。

図４は対象地域の航空写真（オルソ画像）の一例である。図５〜図９は図４に示す対象地域にて取得された航空レーザ計測データ３０に基づく特徴量画像の例であり、図５はＤＣＨＭ画像、図６はＲＩ画像、図７はＢＲＩ画像、図８はＦＰＲ画像、図９はＩＰＲ画像である。

林相区画生成部２４は特徴量画像８０を用いて林相区画生成処理Ｓ１２を行う。既に述べたように、林相区画生成部２４は、特徴量画像生成処理Ｓ１０にて生成された特徴量画像８０上にて特徴量の類似性に基づいて対象地域を複数領域に分割して林相区画を生成する。本実施形態では複数種類の区画用特徴量を用いることに対応して、林相区画生成部２４は区画用特徴量ごとに生成された特徴量画像を複合した複合特徴量画像に対して画像領域分割処理を行う。例えば、複合特徴量画像の画素値は、区画用特徴量ごとに生成された特徴量画像の互いに対応する画素における画素値（特徴量）の組で定義することができる。なお、目的とする林相の区別が１種類の区画用特徴量だけで精度良く行える場合には、当該特徴量の特徴量画像だけを用いて林相区画の生成を行うこともできる。また、複合特徴量画像の画素値として、オルソ画像データを解析して求めたテクスチャ情報やスペクトル情報を併用することも可能である。

初期区画部７０は階層的な領域分割処理の最初の階層の処理を行う。具体的には初期区画部７０は複合特徴量画像（原画像）を複数領域に分割し１次の林相区画を生成する。初期区画部７０は、互いに隣接する複数の画素からなる領域を１次の林相区画とするか否かを、それら画素の画素値の類似性、領域の画像テクスチャの類似性、及び領域の形状や面積などの幾何条件に基づいて決定する。ここで、結合される領域間や結合前後の領域間での画素値の類似性は、画素値の平均値などの統計量などの特徴量に基づいて判断することができる。以下、当該特徴量を、航空レーザ計測データ３０から生成される反射パルス指標等の特徴量と区別するために、画素群色特徴量と称する。

後続区画部７２は階層的な領域分割処理の第２の階層以降の処理を行う。具体的には後続区画部７２は、隣接する林相区画を結合して新たな林相区画を生成する逐次区画処理を１回又は複数回行う。逐次区画処理は、互いに隣接する複数の低次の林相区画について、それらにおける画素値の類似性に関する色条件３４と、それらを結合して得られる高次の林相区画についての形状条件３６と、テクスチャ条件３８と、面積条件４０とに基づいて、それら低次の林相区画を結合するか否かを決める。

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

具体的には、領域併合によるオブジェクトの異質性の変化Δｈは、併合前後における画素群色特徴量の異質性の変化Δｈ_colorと、形状の異質性の変化Δｈ_shapeと、テクスチャの異質性の変化Δｈ_textureと、面積の異質性の変化Δｈ_aeraとから次式によって算出される。

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

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

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

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

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

図１０は林相区画生成部２４の概略の処理フロー図である。林相区画生成処理Ｓ１２は区画用特徴量の特徴量画像８０から生成した複合特徴量画像を領域分割して林相区画（林相区画データ４２）を生成する。上述したように本実施形態では区画用特徴量として反射パルス指標、反射強度及び樹高を利用することに対応して、複合特徴量画像はＦＰＲ（又はＩＰＲ）画像、ＲＩ画像、ＤＣＨＭ画像から構成される。林相区画生成処理Ｓ１２は原画像を領域分割した林相区画を階層的に領域分割して、内部が一様な林相からなる領域に対応した林相区画を生成する。この階層的な領域分割処理は初期区画部７０及び後続区画部７２により行われる。領域分割処理では特徴量画像を構成する画素の画素値を用いた区画処理が少なくとも１回行われる。具体的には、初期区画部７０は原画像を領域分割してスケールが小さい林相区画を生成する（Ｓ１４）。後続区画部７２は林相区画同士を結合することによりスケールが大きくなった林相区画を生成する（Ｓ１６）。

領域分割処理は色条件３４、形状条件３６、テクスチャ条件３８、及び面積条件４０を満たすように行われる。その際、色条件３４と、形状条件３６と、テクスチャ条件３８と、面積条件４０の比重、つまりそれぞれを領域分割に寄与させる度合は調節することができる。

ここで、領域分割処理の階層ごとに、林相区画を生成する際の条件は異なり得る。すなわち、初回の画素を結合して初期の林相区画を生成する際や、低次の林相区画同士を結合して高次の林相区画を生成する際に、スケールパラメータであるΔｈは生成される林相区画のスケールが徐々に大きくなることを可能とするように設定される。なお、本実施形態では階層的な処理を例に説明したが、必ずしも階層的に行わなくてもよい。

図１１は林相区画生成部２４により生成された林相区画データ４２の例を示す画像であり、ＦＰＲ画像、ＲＩ画像及びＤＣＨＭ画像を複合した特徴量画像から生成された林相区画を示している。同図には、図４に示す対象地域にて抽出された林相区画の境界が黒線で描かれている。図３から理解されるように、ＦＰＲには３種類の樹種（ヒノキ林、スギ林及び広葉樹林）間に差が存在し、ＦＰＲのみで樹種の弁別は可能であるが、これにＲＩを組み合わせることで、３つの樹種がより正確に弁別される。これらにＤＣＨＭを加えて林相区画生成部２４による処理を行うことで、樹種の弁別に加え、同じ高さを持つ林相を１つの区画に集約することができる。また、ＤＣＨＭを利用することにより樹齢が異なる同じ樹種の林相を区別できる。よって、ＦＰＲ、ＲＩ及びＤＣＨＭの全部の組み合わせは、それらのうちのいずれか２つの組み合わせ又は１つのみよりも好適な林相区画を与える。

一方、ＦＰＲ及びＲＩを併用しＤＣＨＭを利用しない場合に林相区画生成部２４により生成される林相区画では、同じ高さを持つ林相が１つの区画に集約されることは期待できないものの、樹種の弁別は好適になされる。

また、上述したように、ＦＰＲのみで３種類の樹種の弁別は可能であるので、ＦＰＲのみを利用して林相区画を生成することも可能である。さらに、ＦＰＲとＤＣＨＭとを組み合わせると、同じ高さを持つ林相が１つの区画に集約されるので、ＦＰＲのみの場合より好適な林相区画が得られる。

ＲＩはスギ林を広葉樹林及びヒノキ林から区別できる。よって、ＲＩのみを用いてもスギ林と他の２つの樹種とが区分けされた林相区画を生成することができる。

なお、本実施形態では林相区画生成部２４が特徴量画像を階層的に領域分割して林相区画を生成する例を説明したが、領域分割は特徴量の同一林相での類似性や異なる林相間での差異が区画生成に反映されるものであれば他の手法を用いることもできる。
図２に示すように演算処理装置４は林相区画データ４２の生成を行うと共に、林相特徴量抽出処理Ｓ２０及び林相判別処理Ｓ２２により林相区分図データ４４を生成する。

林相特徴量抽出処理Ｓ２０は、対象地域を航空レーザ計測データ３０及びオルソ画像データ３２に基づき、予め定められた１又は複数種類の林相特徴量を求める。林相特徴量は林相区画ごとに林相を判別することを可能とする特徴量であり、その値や値の分布範囲に林相に応じて差異が生じる量が用いられる。上述した区画用特徴量には航空レーザ計測データ３０から抽出した特徴量だけを用いたが、林相特徴量には航空レーザ計測データ３０から抽出した特徴量と、オルソ画像データ３２から抽出した特徴量との両方を併用する。例えば、航空レーザ計測データ３０から得られる林相特徴量として、反射パルス指標、レーザの反射強度、当該レーザの反射強度を画素値とした対象地域での画像を二値化したパターンから抽出されたテクスチャ特徴量、及び樹高のうちの少なくとも１つを用いることができる。すなわち、レーザ計測データ解析部２０により生成される特徴量画像の画素値が林相特徴量となり、ＤＣＨＭ画像生成部５０、ＲＩ画像生成部５２、ＢＲＩ画像生成部５４及びＰＲ画像生成部５６が生成するＤＣＨＭ画像、ＲＩ画像、ＢＲＩ画像及び、ＦＰＲ画像又はＩＰＲ画像のいずれか、又は全部を林相特徴量に利用することができる。例えば、本実施形態ではＤＣＨＭ画像、ＲＩ画像及び、ＦＰＲ画像又はＩＰＲ画像を林相特徴量として利用する。

またオルソ画像データ３２から得られる林相特徴量として、テクスチャ情報及びスペクトル情報のうちの少なくとも１つを用いることができる。これらテクスチャ情報、スペクトル情報は空中撮影画像解析部２２により生成され、テクスチャ解析部６０が生成するテクスチャ情報及びスペクトル解析部６２が生成するスペクトル情報のいずれか又は全部を用いることができる。

テクスチャ解析部６０は、オルソ画像データ３２からテクスチャ情報として、テクスチャ画像である画像ＢＩ、画像ＢＧＩ及びＬＢＰ画像を生成する。

画像ＢＩはオルソ画像データ３２を二値化処理した二値画像（Binary Image）である。図１２は画像ＢＩの例を示す説明図である。図１２（ａ）〜（ｄ）はそれぞれスギ林、ヒノキ林、広葉樹林、非森林の例であり、左右に並ぶ２つの画像のうち左側がオルソ画像データ３２、右側が画像ＢＩである。当該画像ＢＩにおいて白領域がオルソ画像データ３２にて輝度がしきい値以上の領域であり、黒領域がしきい値未満の領域である。ちなみに、ここでの二値化のしきい値は例えば大津の手法により決定することができる。例えば、画像ＢＩには日向及び日陰の分布パターンが現れる。

画像ＢＧＩはオルソ画像データ３２における輝度勾配を二値化処理した画像（二値化勾配画像： Binary Gradient Image）である。図１３は画像ＢＧＩの例を示す説明図である。図１３（ａ）〜（ｄ）は図１２と同様、それぞれスギ林、ヒノキ林、広葉樹林、非森林の例であり、左右に並ぶ２つの画像のうち左側がオルソ画像データ３２、右側が画像ＢＧＩである。当該画像ＢＧＩにおいて白領域が輝度勾配がしきい値以上の領域であり、黒領域がしきい値未満の領域である。ここでも、しきい値は例えば大津の手法により決定することができる。例えば、画像ＢＧＩからは日向と日陰とが切り替わる空間的な頻度の多寡が読み取れる。

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

スペクトル解析部６２はオルソ画像データ３２からスペクトル情報を抽出する。本実施形態ではスペクトル解析部６２はＲ，Ｇ，Ｂ，ＮＩＲの４バンドのマルチスペクトル画像に対して次式で示される正規化処理を行い、正規化後の成分Ｒ’，Ｇ’，Ｂ’からなる画像を生成する。

なお、スペクトル解析部６２は正規化処理を行わないＲ，Ｇ，Ｂ，ＮＩＲの４バンドのマルチスペクトル画像を出力する構成とすることもできる。

上述のようにＤＣＨＭ画像生成部５０、ＲＩ画像生成部５２、ＢＲＩ画像生成部５４、ＰＲ画像生成部５６、テクスチャ解析部６０、スペクトル解析部６２が林相特徴量抽出処理Ｓ２０を行う林相特徴量抽出手段として機能し、林相特徴量データ８２を生成する。ここで、レーザ計測データ解析部２０を構成する複数の画像生成部の中には上述したように、区画用特徴量に関する特徴量画像８０を生成する特徴量画像生成手段を兼ねるものがあり、当該画像生成部が生成する特徴量画像は、特徴量画像生成処理Ｓ１０及び林相特徴量抽出処理Ｓ２０それぞれの処理結果として共用することができる。そこで、本実施形態では、特徴量画像８０のうち共用可能なものを林相特徴量抽出処理Ｓ２０が林相特徴量データ８２として利用する構成としている。

林相区分決定部２６は既に述べたように林相判別手段であり、林相特徴量データ８２及び林相区画データ４２を用い、林相区画ごとの林相特徴量に基づいて、当該区画における林相（林相区分）を判別する林相判別処理Ｓ２２を行う。例えば、林相区分決定部２６は最近隣法（Nearest Neighbor Method）を用いて林相を判別する。

最近隣法は教師付き分類の１種であり、教師付き分類は教師データ（トレーニングデータ）と呼ばれる予め分類結果が定まっているデータを用いて未分類データを分類する。例えば、分類対象のオブジェクトは、教師データに基づいて定義されるクラスのうち、当該オブジェクトに最も近いオブジェクトを有するクラスに割り当てられる。

林相区分決定部２６は、複数種類の林相特徴量それぞれについて各林相区画内における値（画素値）を代表する量（区画毎特徴量）を算出し、林相区画ごとに複数種類の区画毎特徴量の組で表される座標値を定義する。すなわち、区画毎特徴量の種類の数をＭとすると、Ｍ次元空間内にて、分類対象となるオブジェクトである林相区画の座標が定義される。例えば、或る林相特徴量についての区画毎特徴量は、林相区画内での当該林相特徴量の平均値、標準偏差などの統計量で定義することができる。最近隣法におけるオブジェクト間の距離は当該Ｍ次元空間での距離で定義することができる。

教師データはユーザにより与えられる。ユーザは例えば、林相区画データ４２で与えられる林相区画のうち林相を判読したものを教師データとして選択する。林相区分決定部２６は教師データを用いて、林相区画の林相を最近隣法により判別し林相区分図（林相区分図データ４４）を作成する。

図１５は図４に示す対象地域での最近隣法による林相分類で作成された林相区分図であり、ヒノキ林区画９０、スギ林区画９２、広葉樹林区画９４及び非森林区画９６に分類されている。

ここでは、林相区画の分類手法として最近隣法を説明したが、これは一例であり、航空レーザ計測データ３０及びオルソ画像データ３２から得られる林相特徴量に基づく林相判別をテンプレートマッチング法（Template Matching Method）や決定木法（Decision Tree Method）など、最近隣法以外の手法で行うこともできる。例えば、テンプレートマッチング法はテンプレートと呼ばれる予め分類結果が定まっているデータを用いて未分類データを分類する。例えば、対象地域の地点（画素）ごとに複数種類の林相特徴量それぞれを座標軸とした多次元空間での座標値を画素値として定義し、林相区画内の画素値のヒストグラムについてヒストグラム間の距離に基づいてマッチングを行い、分類対象を距離が最も近いテンプレートの林相に分類する。ヒストグラムの距離は例えば、Matsusita Distanceを用いて測ることができる。

以上説明したように、本発明による林相解析システム２は航空写真等の空中撮影画像とレーザ計測データとを併用して林相区分図を自動生成する。

例えば、航空写真からは、広い波長帯域の色情報が対象地域にて高解像度に得られる一方、色は撮影季節の影響を受けたり、異なる樹種が色的に非常に類似する場合が生じたりする。一方、レーザ計測データからは森林の地形や樹高、樹冠表面及び内部の構造などに関する情報が得られる一方、航空写真ほどには高解像度の情報は得られず、また近赤外波長領域での情報しか得られない。林相解析システム２ではこれら航空写真とレーザ計測データとを併用することで、それぞれの特長を活かして高精度の林相区分図を自動生成することが可能となる。

例えば、レーザ計測データから得られる高さ情報を用いることにより、同じ樹種でも樹齢や生育状況が異なる林相を区別することができる。一方、航空写真からは高解像度のテクスチャ情報が生成され林相の区別に活用される。また、レーザ計測データを併用することで、撮影の時期、撮影の方向・角度、及び太陽の方位角・高度角などにより航空写真が受ける影響が林相解析に与える影響が相対的に軽減され、より安定かつ正確に林相区分図を作成することが可能となる。また、レーザ計測データを併用することで、航空写真より低分解能である衛星画像を用いて林相区分図を作成することが容易となる。例えば、レーザ計測データから抽出したＲＩ画像を二値化したＢＲＩ画像をテクスチャ情報として用いることで、衛星画像から抽出されるテクスチャ情報が低分解能となることを補うことができる。衛星画像を利用することで撮影コストの軽減が可能となり、林相区分図の作成・更新をより効率的に行うことができる。

２林相解析システム、４演算処理装置、６記憶装置、８入力装置、１０出力装置、２０レーザ計測データ解析部、２２空中撮影画像解析部、２４林相区画生成部、２６林相区分決定部、３０航空レーザ計測データ、３２オルソ画像データ、３４色条件、３６形状条件、３８テクスチャ条件、４０面積条件、４２林相区画データ、４４林相区分図データ、５０ＤＣＨＭ画像生成部、５２ＲＩ画像生成部、５４ＢＲＩ画像生成部、５６ＰＲ画像生成部、６０テクスチャ解析部、６２スペクトル解析部、７０初期区画部、７２後続区画部、８０特徴量画像、８２林相特徴量データ。

Claims

森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成手段と、
前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成手段と、
前記空中撮影画像に基づいて少なくとも１種類の林相特徴量を求め、かつ前記航空レーザ計測データに基づいて少なくとも１種類の林相特徴量を求める林相特徴量抽出手段と、
前記林相区画ごとの複数種類の前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別手段と、
を有することを特徴とする林相解析装置。
請求項１に記載の林相解析装置において、
前記区画用特徴量は、レーザの照射パルスに対する反射パルスの個数に応じて定まる前記各計測地点での反射パルス指標、レーザの反射強度及び樹高のうちの少なくとも１つを含むこと、を特徴とする林相解析装置。
請求項１又は請求項２に記載の林相解析装置において、
前記林相特徴量は、レーザの照射パルスに対する反射パルスの個数に応じて定まる前記各計測地点での反射パルス指標、レーザの反射強度、当該レーザの反射強度を画素値とした前記対象地域での画像を二値化したパターンから抽出されたテクスチャ特徴量、及び樹高のうちの少なくとも１つと、それぞれ前記空中撮影画像から抽出されるテクスチャ情報及びスペクトル情報のうちの少なくとも１つとを含むこと、
を特徴とする林相解析装置。
森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成ステップと、
前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成ステップと、
前記空中撮影画像に基づいて少なくとも１種類の林相特徴量を求め、かつ前記航空レーザ計測データに基づいて少なくとも１種類の林相特徴量を求める林相特徴量抽出ステップと、
前記林相区画ごとの複数種類の前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別ステップと、
を有することを特徴とする林相解析方法。
コンピュータに林相解析を行わせるためのプログラムであって、当該コンピュータを、
森林を含む対象地域をスキャンして取得した航空レーザ計測データ、又は当該航空レーザ計測データ及び前記対象地域を上空から撮影した空中撮影画像に基づき、予め定められた１又は複数種類の区画用特徴量を各計測地点について求め、当該区画用特徴量を画素値とする前記対象地域の特徴量画像を生成する特徴量画像生成手段、
前記区画用特徴量の類似性に基づいて前記特徴量画像上にて前記対象地域を複数領域に分割して林相区画を生成する林相区画生成手段、
前記空中撮影画像に基づいて少なくとも１種類の林相特徴量を求め、かつ前記航空レーザ計測データに基づいて少なくとも１種類の林相特徴量を求める林相特徴量抽出手段、及び、
前記林相区画ごとの複数種類の前記林相特徴量に基づいて、当該林相区画における林相を判別する林相判別手段、
として機能させることを特徴とするプログラム。