葉面積指数算出装置、葉面積指数算出方法及びそのプログラム
【課題】航空レーザー計測を用いて葉面積指数を高精度に推定することができる葉面積指数算出方法を提供する。
【解決手段】上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手順と、地盤の反射率を記憶する手順と、記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手順（Ｓ１３）と、記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手順（Ｓ１４）と、記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出する手順（Ｓ１２）と、樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、及び樹冠で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する手順（Ｓ１５）とを含む。
【解決手段】上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手順と、地盤の反射率を記憶する手順と、記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手順（Ｓ１３）と、記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手順（Ｓ１４）と、記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出する手順（Ｓ１２）と、樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、及び樹冠で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する手順（Ｓ１５）とを含む。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、林分パラメータとして重要な葉面積指数を算出する葉面積指数算出装置、葉面積指数算出方法及びそのプログラムに関する。
【背景技術】
【０００２】
樹木本数、樹冠面積、樹高、胸高直径、及び葉面積指数(Leaf Area Index:ＬＡＩ)等の林分パラメータを把握することは、森林の状態を定量化し、林業や保全活動へ反映させる上で重要である。ここで、ＬＡＩは、植物の葉の多少の度合いを示す指標であり、地上の単位面積に対しての、その上方に存在するすべての葉の片側の総面積の比率である（例えば、特許文献１参照。）。
【０００３】
ＬＡＩ推定手法としては、現地調査を行うのが一般的である。しかし、この方法では、局所的にしか推定することができず、広域のＬＡＩを推定するのには、時間や労力がかかる。
【０００４】
そこで、広域にＬＡＩを推定する場合には、衛星データを活用する手法がある。これは、衛星データより正規化植生指標（ＮＤＶＩ）を求め、ＮＤＶＩとＬＡＩの回帰分析を実施することによって得られた回帰式に基づいて推定をするものである。この手法においても、ＬＡＩの実測値を取得していることが前提であり、広域展開できるものの、実測データがなければ実施できない。
【０００５】
また、航空レーザーを活用して広域のＬＡＩを推定する手法がある。この方法では、航空レーザーにより計測されたレーザー反射地点の高さとＤＥＭの差をとった地上面からの高さ（樹高データ）を用いてＬＡＩを推定する。具体的には、横軸を地上面からの高さ、縦軸を頻度としたヒストグラムから累積頻度分布を作成し、これが光の減衰を表現していると仮定してＬＡＩを算出している。
【特許文献１】特開２００５−５２０４５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来の航空レーザーを活用したＬＡＩ推定手法では、ヒストグラムが樹冠内部の情報を反映しきれていないため、算出したＬＡＩには疑問が残る。即ち、ＬＡＩを高精度に推定することができない問題点があった。
【０００７】
上記問題点を鑑み、本発明は、航空レーザー計測を用いて広域かつ簡易にＬＡＩを高精度に推定することができる葉面積指数算出装置、葉面積指数算出方法及びそのプログラムを提供することである。
【課題を解決するための手段】
【０００８】
本願発明の一態様によれば、（イ）上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶するステップと、（ロ）地盤の反射率を記憶するステップと、（ハ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出するステップと、（ニ）記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出するステップと、（ホ）記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出するステップと、（ヘ）記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出するステップと、（ト）樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出するステップとを含む葉面積指数算出方法が提供される。
【０００９】
本願発明の他の態様によれば、コンピュータを、（イ）上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手段と、（ロ）地盤の反射率を記憶する手段と、（ハ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手段と、（ニ）記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手段と、（ホ）記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出する手段と、（ヘ）記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出する手順と、（ト）樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する手段手段として機能させるためのプログラムが提供される。
【００１０】
本願発明の更に他の態様によれば、（イ）上空から地上へレーザーを照射するレーザー照射手段と、（ロ）地上で反射したレーザーの反射強度データを記憶する反射強度記憶手段と、（ハ）地盤の反射率を記憶する反射率記憶手段と、（ニ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出し、記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出し、記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出し、記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出し、樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する演算手段と、（ホ）葉面積指数を出力する出力手段とを備える葉面積算出装置が提供される。
【発明の効果】
【００１１】
本発明によれば、航空レーザー計測を用いてＬＡＩを高精度に推定することができる葉面積指数算出装置、葉面積指数算出方法及びそのプログラムを提供することができる。
【００１２】
したがって、このＬＡＩを林分パラメータとして用いて、ＣＯ2削減量を精度良く推定すること等が可能となる。
【発明を実施するための最良の形態】
【００１３】
次に、図面を参照して、本発明の第１及び第２の実施の形態を説明する。以下の図面の記載において、同一又は類似の部分には同一又は類似の符号を付している。また、以下に示す実施の形態は、この発明の技術的思想を具体化するための装置や方法を例示するものであって、この発明の技術的思想は、構成部品の材質、形状、構造、配置等を下記のものに特定するものでない。この発明の技術的思想は、特許請求の範囲において、種々の変更を加えることができる。
【００１４】
本発明の第１の実施の形態に係る葉面積指数算出装置は、図１に示すように、中央処理装置（ＣＰＵ）１、表示部２、入力部３、出力部４及び記憶装置５を備える。入力部３としては、例えばキーボード、マウス、ＯＣＲ等の認識装置、スキャナ、カメラ等の画像入力装置、マイク等の音声入力装置等が使用可能である。出力部４としては、液晶ディスプレイ、ＣＲＴディスプレイ等の表示装置や、インクジェットプリンタ、レーザープリンタ等の印刷装置等を用いることができる。表示部２は、液晶ディスプレイ、ＣＲＴディスプレイ等である。記憶装置５には、ＲＯＭ及びＲＡＭが組み込まれている。ＲＯＭは、ＣＰＵ１において実行されるプログラムを格納する。ＲＡＭは、ＣＰＵ１におけるプログラム実行処理中に利用されるデータ等を一時的に格納したり、作業領域として利用される一時的なデータメモリ等として機能する。記憶装置５としては、例えば半導体メモリ、磁気ディスク、光ディスク、光磁気ディスクや磁気テープ等が採用可能である。
【００１５】
図２は、本発明の第１の実施の形態の葉面積指数算出装置の概略構成図である。（第１の）地盤反射強度抽出部６、入射強度算出部７、表示処理部１０、（第２の）地盤反射強度抽出部１１、反射率算出部１２、樹冠反射強度抽出部１３、樹冠通過強度抽出部１４、及び葉面積指数算出部１５は、図１に示したＣＰＵ１により実現される。
【００１６】
図２に示す地盤反射強度記憶部８、反射率記憶部９、オルソフォト画像記憶部２０、反射強度データ記憶部２１、反射率記憶部２２、樹木頂点記憶部２３、ＤＥＭ記憶部２４、閾値記憶部２５、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８、樹冠通過強度記憶部２９及び葉面積指数記憶部３１は、図１に示した記憶装置５に含まれる。
【００１７】
オルソフォト画像記憶部２０は、ＬＡＩ算出対象となる森林域Ａを含むオルソフォト画像を格納する。「オルソフォト画像」とは、写真画像に三次元計測データ等を与えて正射変換を行った画像である。オルソフォト画像は、例えば、座標（Ｘ，Ｙ）とそれに対応する赤（Ｒ）、緑（Ｇ）、青（Ｂ）の色情報を含む。
【００１８】
反射強度データ記憶部２１は、航空レーザー（以下、単に「レーザー」という。）の反射強度データを格納する。反射強度データは、航空機に設けたレーザー照射手段（図示省略）によりレーザーを地上へ照射し、航空機内のレーザー受信手段（図示省略）により地上で反射したレーザーを受信することで取得することができる。反射強度データは、任意の地点で得られるランダムデータであり、３次元座標（Ｘ，Ｙ，Ｚｄ）とレーザーの反射強度を含む。高さの座標Ｚｄは、地盤や樹木等のレーザーの反射地点の高さを示す。
【００１９】
レーザーの反射強度は、パルス波形で示される。パルス波形には複数のピーク値があり、最初のピーク値を「ファーストパルスの強度Ifirst」、最後のピーク値を「ラストパルスの強度Ilast」と定義する。
【００２０】
ＤＥＭ記憶部２４は、森林の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を格納する。ＤＥＭは、例えば平面直角座標系のＸ軸、Ｙ軸を一定間隔（０．５ｍ）で区切ったメッシュの格子点に割り当てている。ＤＥＭは、３次元座標（Ｘ，Ｙ，Ｚh）を含む。高さの座標Ｚｈは、地盤の高さを示す。
【００２１】
地盤反射強度抽出部６は、反射強度データ記憶部２１から、既知反射率領域内の座標（Ｘ，Ｙ）の複数の反射強度データに含まれるファーストパルスの強度Ifirstをそれぞれ読み出す。ここで、「既知反射率領域」とは、地盤の反射率が予め分かっている領域をいい、例えばアスファルトや草原の領域である。更に、地盤反射強度抽出部６は、複数のファーストパルスの強度Ifirstの平均値を算出し、その平均値を、図３に示すような、既知反射率領域における地盤に直接達して反射したレーザーの反射強度（以下、「地盤反射強度」という。）Iaとして抽出する。既知反射率領域における地盤反射強度Iaは、地盤反射強度記憶部８に記憶される。
【００２２】
入射強度算出部７は、地盤反射強度記憶部８から、既知反射率領域における地盤反射強度Iaを読み出す。更に、入射強度算出部７は、反射率記憶部９から反射強度Iaに対応する既知の反射率Ｒaを読み出す。既知の反射率Ｒaは、図４に示すような一般的な分光反射特性のグラフから設定可能である。入射強度算出部７は、既知反射率領域における地盤反射強度Ia及び既知の反射率Ｒaを用いて、式（１）のように、図３に示すような、地盤に入射するレーザーの強度（以下、「入射強度」という。）I0を算出する。
【００２３】
I0＝Ia／Ra …（１）
例えば、地盤反射強度Iaの値を１２．２とし、使用したレーザーの波長が１．０６４［μm］程度であることを考慮して反射率Ｒaの値を０．１４とすると、入射強度I0の値は、８７．１２となる。算出された入射強度I0は、入射強度記憶部２７に格納される。
【００２４】
樹木頂点記憶部２３は、森林域Ａ内の樹木毎に、樹木頂点の座標を格納する（この樹木頂点の算出方法については後述する）。本発明の第１の実施の形態においては、一例として、図５に示すオルソフォト画像において、森林域Ａ内の樹木（樹冠）Ｔ１〜Ｔ６のＬＡＩを算出することとする。例えば、オルソフォト画像記憶部２０が、樹木（樹冠）Ｔ１〜Ｔ６の形状・座標等を、オルソフォト画像に対応させて格納していれば良い。なお、ＬＡＩの算出対象は樹木Ｔ１〜Ｔ６に限定されず、例えば森林域Ａ内のすべての樹木のＬＡＩを算出しても良い。樹木頂点記憶部２３は、図５に示した樹木Ｔ１〜Ｔ６毎に、樹木頂点Ｐi（ｉ＝１〜６）の座標（Ｘi，Ｙi）を格納する。
【００２５】
樹冠反射強度抽出部１３は、樹木頂点記憶部２３から樹木頂点Ｐiの座標（Ｘi，Ｙi）を読み出し、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）のラストパルスの強度Ifirstを、図３に示すような樹冠で反射したレーザーの反射強度（以下、「樹冠反射強度」という。）Ictとして抽出する。樹冠反射強度Ictは、樹冠反射強度記憶部２８に格納される。
【００２６】
樹冠通過強度抽出部１４は、樹木頂点記憶部２３に格納された樹木頂点Ｐiの座標（Ｘi，Ｙi）を参照して、反射強度データ記憶部２１から、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データの反射地点の高さ（Ｚ値）Ｚdを読み出す。更に、樹冠通過強度抽出部１４は、ＤＥＭ記憶部２４から、ＤＥＭの樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の地盤の高さＺhを読み出す。更に、樹冠通過強度抽出部１４は、反射地点の高さＺdと地盤の高さＺhとの差分値Ｚkを算出する。
【００２７】
更に、樹冠通過強度抽出部１４は、閾値記憶部２５から閾値（例えば、１ｍ）を読み出して、差分値Ｚkと比較する。差分値Ｚkが閾値以内の場合、レーザーが樹冠を通過し、地盤で反射したものと判断されるので、反射強度データのラストパルスの強度Ilastを、図３に示すような樹冠を通過して減衰した後に地盤で反射したレーザーの反射強度（以下、「樹冠通過強度」という。）Icbとして抽出する。図３に示すように、樹冠通過強度Icbは、樹冠を通過することで減衰しているので、樹冠内部の情報を有するといえる。抽出した樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。
【００２８】
地盤反射強度抽出部１１は、反射強度データ記憶部２１から、森林域Ａ内の反射強度データの反射地点の高さ（Ｚ値）Ｚdを読み出す。更に、地盤反射強度抽出部１１は、ＤＥＭ記憶部２４から、森林域Ａ内のＤＥＭの地盤の高さＺhを読み出す。更に、地盤反射強度抽出部１１は、反射地点の高さＺdと地盤の高さＺhとの差分値Ｚkを算出する。
【００２９】
更に、地盤反射強度抽出部１１は、閾値記憶部２５から閾値（例えば、１ｍ）を読み出して、差分値Ｚkと比較する。差分値Ｚkが閾値よりも小さい場合、レーザーが直接地盤に到達したものと判断されるので、反射強度データのファーストパルスの強度Ifirstを、図３に示すような地盤に直接達して反射したレーザーの反射強度（以下、「地盤反射強度」という。）Ibとして抽出する。抽出した樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。なお、樹冠通過強度Icbが森林域Ａ内の反射率Rbで得られることから、地盤反射強度Ibも森林域Ａから抽出している。
【００３０】
反射率算出部１２は、入射強度記憶部２７から入射強度I0を読み込み、地盤反射強度記憶部２６から森林域Ａ内の地盤反射強度Ibを読み込んで、式（２）を用いて、森林域Ａ内の反射率Ｒｂを算出する。
【００３１】
Rb＝Ib／I0 ・・・（２）
森林域Ａ内の反射率Ｒｂは、反射率記憶部２２に記憶される。
【００３２】
図６に示すように、樹木Ｔ１〜Ｔ６毎に、森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0が設定される。
【００３３】
葉面積指数算出部１５は、地盤反射強度記憶部２６、入射強度記憶部２７、及び樹冠反射強度記憶部２８、樹冠通過強度記憶部２９のそれぞれに格納された森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出して、樹木Ｔ１〜Ｔ６毎に、ＬＡＩを算出する。
【００３４】
地盤反射強度Ibは下の式で表される。
【００３５】
Ib＝I0・Rｂ …（３）
ベール−ランベルトの法則を参考にすると、葉間隙（gap fraction）は下の式（４）で表される。
【００３６】
Ｔ(ｚ)＝exp(−Ｇ・μ・ｚ) …（４）
ここで、Ｔは葉間隙、Ｇは葉の投影係数、μは葉群密度、ｚは梢端からの透過距離を表す。
【００３７】
樹冠通過強度Icbは、下の式（５）で表される。
【００３８】
Icb＝(I0−Ict) ・exp(-2G・μ・d)・Rb …（５）
ここで、ｄは樹冠厚さを表す。式（５）より、樹冠通過強度Icbと 森林域Ａ内の地盤反射強度Ibの比をとると式（６）のようになる。
【００３９】
Icb／Ib＝ ((I0−Ict)/ I0) ・ exp(-2G・μ・d) …（６）
葉面積指数は下の式（７）で表される。
【００４０】
ＬＡＩ＝μ・ｄ …（７）
葉を水平面に投影したと仮定して、Ｇ＝１とする。式（６）及び式（７）より、葉面積指数ＬＡＩは、式（８）で求めることができる。
【００４１】
ＬＡＩ＝-1/2 ・ ln( (Icb・I0) / (Ib・(I0 - Ict)) ) …（８）
葉面積指数算出部１５は、樹木Ｔ１〜Ｔ６毎に、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８及び樹冠通過強度記憶部２９から、森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出す。更に、葉面積指数算出部１５は、式（８）のように、ＬＡＩを算出する。算出したＬＡＩは、葉面積指数記憶部３１に格納される。葉面積指数記憶部３１は、図７に示すように、樹木Ｔ１〜Ｔ６毎に、ＬＡＩの値を格納する。出力部４は、ＬＡＩの値を、記録媒体に記録したり、画面に表示したり、印刷したりする。
【００４２】
次に、本発明の第１の実施の形態に係る葉面積指数算出方法の一例を、図８のフローチャートを参照しながら説明する。
【００４３】
（イ）ステップＳ１１において、地盤反射強度抽出部１１が、反射強度データ記憶部２１から、図３に示すような地盤に直接達して反射したレーザーの反射強度（地盤反射強度）Ibを抽出する。森林域Ａ内の地盤反射強度Ibは、地盤反射強度記憶部２６に格納される。
【００４４】
（ロ）ステップＳ１２において、地盤反射強度抽出部６が、反射強度データ記憶部２１から、図３に示すような既知反射率領域の地盤反射強度Iaを抽出し、地盤反射強度記憶部８に記憶させる。そして、入射強度算出部７が、地盤反射強度記憶部８から、既知反射率領域の地盤反射強度Iaを読み出し、更に反射率記憶部９から既知の反射率Ｒaを読み出して、式（１）を用いて、樹木Ｔ１〜Ｔ６の樹冠に到達する直前のレーザーの入射強度（入射強度）I0を算出する。入射強度I0は、入射強度記憶部２７に格納される。
【００４５】
（ハ）ステップＳ１３において、樹冠反射強度抽出部１３が、反射強度データ記憶部２１から、図３に示すような樹木Ｔ１〜Ｔ６の樹冠で反射したレーザーの反射強度（樹冠反射強度）Ictを、樹木Ｔ１〜Ｔ６毎に抽出する。樹冠反射強度Ictは、樹冠反射強度記憶部２８に格納される。
【００４６】
（ニ）ステップＳ１４において、樹冠通過強度抽出部１４は、反射強度データ記憶部２１から、図３に示すような、樹木Ｔ１〜Ｔ６の樹冠を通過した後に地盤で反射したレーザーの反射強度（樹冠通過強度）Icbを、樹木Ｔ１〜Ｔ６毎に抽出する。樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。
【００４７】
（ホ）ステップＳ１５において、葉面積指数算出部１５が、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８及び樹冠通過強度記憶部２９から森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出して、式（７）を用いて、樹木Ｔ１〜Ｔ６毎にＬＡＩを算出する。ＬＡＩは、葉面積指数記憶部３１に格納される。
【００４８】
なお、ステップＳ１１〜Ｓ１４の手順は、それぞれ独立して行うことが可能であり、図８のフローチャートに示した順序に特に限定されない。
【００４９】
次に、図８のステップＳ１１の森林域Ａ内の地盤反射強度Ibの設定方法の詳細を、図９のフローチャートを参照しながら説明する。
【００５０】
（イ）ステップＳ２１において、表示部２が、オルソフォト画像を画面に表示している。利用者が、オルソフォト画像の森林域ＡをＬＡＩの算出対象とする指示を入力部３に入力する。表示処理部１０が、入力部３からの指示に応じて、森林域Ａを示す画面上のピクセル座標（Ｘｐ，Ｙｐ）を読み出す。そして、表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当するオルソフォト画像上の平面直角座標（Ｘo，Ｙo）を読み、地盤反射強度抽出部１１に引き渡す。
【００５１】
（ロ）ステップＳ２２において、地盤反射強度抽出部１１が、森林域Ａ内の反射強度データを引き当てる。ステップＳ２３において、地盤反射強度抽出部１１が、反射強度データに含まれるレーザーの反射地点の高さＺｄを抽出する。
【００５２】
（ハ）ステップＳ２４において、地盤反射強度抽出部１１が、ＤＥＭ記憶部２４から、森林域Ａ内のＤＥＭの地盤の高さＺｈを読み出す。ステップＳ２５において、地盤反射強度抽出部１１が、反射強度データの反射地点の高さＺｄと、ＤＥＭの地盤の高さＺｈとの差分値Ｚｋをとる。
【００５３】
（ニ）ステップＳ２６において、地盤反射強度抽出部１１が、差分値Ｚｋが閾値と比較する。差分値Ｚｋが閾値以上の場合は、ステップＳ２４に戻り、森林域Ａ内の他の反射強度データを引き当てる。他方、差分値Ｚｋが閾値以内場合は、ステップＳ２７に進む。
【００５４】
（ホ）ステップＳ２７において、地盤反射強度抽出部１１が、反射強度データのファーストパルスの強度Ifirstを、図３に示すような地盤反射強度Ibとして抽出する。ステップＳ２８において、地盤反射強度抽出部１１が、地盤反射強度Ibを地盤反射強度記憶部２６に記憶させる。
【００５５】
次に、図８のステップＳ１３の樹冠反射強度Ictの設定方法の詳細を、図１０のフローチャートを参照しながら説明する。
【００５６】
（イ）ステップＳ３１において、表示部２が、オルソフォト画像を表示している。表示処理部１０が、入力部３からの指示に応じて、表示部２の画面に表示された森林域Ａを示すピクセル座標（Ｘｐ，Ｙｐ）を読み出す。表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当するオルソフォト画像の平面直角座標（Ｘo，Ｙo）を読み、樹冠反射強度抽出部１３に引き渡す。
【００５７】
（ロ）ステップＳ３２において、樹冠反射強度抽出部１３が、平面直角座標（Ｘo，Ｙo）に基づいて、樹木頂点記憶部２３から森林域Ａ内に存在する樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）を読み出す。
【００５８】
（ハ）ステップＳ３３において、樹冠反射強度抽出部１３が、樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）のうち、樹木頂点Ｐ１に着目する。
【００５９】
（ニ）ステップＳ３４において、樹冠反射強度抽出部１３が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致、或いは最も近い座標（Ｘ，Ｙ）を有する反射強度データのファーストパルスの強度Ifirstを樹冠反射強度Ictとして抽出する。
【００６０】
（ホ）ステップＳ３５において、樹冠反射強度抽出部１３が、抽出した樹冠反射強度Ictを、樹木Ｔ１の樹木頂点Ｐ１に対応させて樹冠反射強度記憶部２８に記憶させる。
【００６１】
（ヘ）ステップＳ３６において、樹冠反射強度抽出部１３が、着目した樹木頂点Ｐｉが最後の樹木頂点Ｐ６であるか、即ち全ての樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictを設定したか判断する。全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictを抽出していなければ、ステップＳ３７においてｉを１インクリメントし、ステップＳ３４に戻り、次の樹木頂点Ｐ２に対して樹冠反射強度Ictを抽出する。この結果、全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictが抽出される。
【００６２】
次に、図８のステップＳ１４の樹冠通過強度Icbの設定方法の詳細を、図１１のフローチャートを参照しながら説明する。
【００６３】
（イ）ステップＳ４１において、表示部２が、オルソフォト画像を表示している。表示処理部１０が、入力部３からの指示に応じて、表示部２の画面に表示されたＬＡＩの算出対象となる樹木群の森林域Ａを示すピクセル座標（Ｘｐ，Ｙｐ）を読み出す。表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当する平面直角座標（Ｘo，Ｙo）を読み、樹冠通過強度抽出部１４に引き渡す。
【００６４】
（ロ）ステップＳ４２において、樹冠通過強度抽出部１４が、樹木頂点記憶部２３から、森林域Ａ内に存在する樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）を読み出す。
【００６５】
（ハ）ステップＳ４３において、樹冠通過強度抽出部１４が、読み出した樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）のうち、樹木頂点Ｐ１に着目する。
【００６６】
（ニ）ステップＳ４４において、樹冠通過強度抽出部１４が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致、或いは最も近い座標（Ｘ，Ｙ）を有する反射強度データを引き当てる。
【００６７】
（ホ）ステップＳ４５において、樹冠通過強度抽出部１４が、反射強度データの反射地点の高さＺｄを抽出する。
【００６８】
（ヘ）ステップＳ４６において、樹冠通過強度抽出部１４が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致する、或いは最も近いＤＥＭの地盤の高さＺｈを読み出す。
【００６９】
（ト）ステップＳ４７において、樹冠通過強度抽出部１４が、反射強度データの反射地点の高さＺｄと、ＤＥＭの地盤の高さＺｈとの差分値Ｚｋをとる。
【００７０】
（チ）ステップＳ４８において、樹冠通過強度抽出部１４が、差分値Ｚｋが、閾値以上か判断する。閾値よりも大きい場合は、ステップＳ４９に進み樹木頂点Ｐｉの座標（Ｘｉ，Ｙｉ）に次に近い座標を有する反射強度データを引き当て、ステップS４４に戻る。他方、閾値以内の場合は、ステップＳ５０に進み、樹冠通過強度抽出部１４が、Ｚｄを有する反射強度データを引き当てる。
【００７１】
（リ）ステップＳ５１において、樹冠通過強度抽出部１４が、反射強度データのラストパルスの強度Ilastを樹冠通過強度Icbとして抽出する。
【００７２】
（ヌ）ステップＳ５２において、樹冠通過強度抽出部１４が、樹冠通過強度Icbを樹冠通過強度記憶部２９に樹木頂点Ｐ１と対応させて保存（設定）する。
【００７３】
（ル）ステップＳ５３において、樹冠通過強度抽出部１４が、着目したのが樹木頂点Ｐ６であるか、即ち全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbを設定したか判断する。全ての樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbを設定していなければ、ステップＳ５４に進み、次の樹木頂点Ｐ２を対象として、ステップＳ４４に戻る。この結果、全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbが抽出される。
【００７４】
本発明の第１の実施の形態によれば、従来、現地調査により局所的なＬＡＩを求めることしかできなかったが、航空レーザー計測でＬＡＩを推定できるので、現地に立ち入ることなく、広域かつ迅速にＬＡＩを推定することができる。
【００７５】
更に、地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0等のレーザーの強度を活用し、樹冠内部の情報（レーザー反射強度の減衰量）を用いてＬＡＩを推定しているので、高精度なＬＡＩを推定することができる。
【００７６】
更に、樹木頂点記憶部２３に格納された樹木頂点Ｐiを用いて、樹冠反射強度抽出部１３が、樹木頂点Ｐiに一致、或いは最も近い座標（Ｘ，Ｙ）の樹冠反射強度Ictを抽出し、樹冠通過強度抽出部１４が、樹木頂点Ｐiに一致、或いは最も近い座標（Ｘ，Ｙ）の樹冠通過強度Icbを抽出することにより、ＬＡＩを推定することができる。樹木頂点Ｐiの決定方法は後述する。
【００７７】
このようにして求めたＬＡＩは森林のＣＯ2吸収量を推定する際のパラメータとして活用することもできる。
【００７８】
（第２の実施の形態）
本発明の第２の実施の形態として、樹木頂点の特定方法について説明する。図１３に示した樹木頂点記憶部２３に含まれる樹木頂点Ｐｉは、以下のように特定している。
【００７９】
本発明の第２の実施の形態は、樹冠を際立たせて再現し、単木の頂点を認識するものである。この単木の頂点の認識のために、新たに「尾根谷度」データを作成することによって樹冠形状を判別する。
【００８０】
尾根谷度は次の様に算出する。
【００８１】
本発明の第２の実施の形態では、ＤＥＭとＤＳＭの差のデータである樹木高のグリッドデータ（0.5m）を用いて、８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。検索範囲は1m・2m・5m・10mの4通り試算し、樹冠形状を適切に表現できた2mに決定した。そして、検索範囲2mにおける大まかな地形を表す、地上開度と地下開度とに挟まれる角度の2等分線とを尾根谷度（φ3）と定義する。
【００８２】
以下に図面を用いて樹木認識の実施の形態を説明する。
【００８３】
図１２は本発明の第２の実施の形態の樹木頂点認識装置の概略構成図である。図１２に示すように、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を記憶したデータベース１１０と、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model）を記憶したデータベース１１１と、特定地域のデジタル写真を記憶したデータベース１１２と、表示部１３と、樹木高算出部１１４と、樹木高データが記憶されるファイル（メモリ）１１５と、平滑フィルタ部１１６と（３ピクセル×３ピクセル）、尾根谷度算出部１１７と、尾根谷度が記憶されるファイル１１８と、表示処理部１１９と、局所最大フィルタ部１２０と、樹木頂点候補抽出部１２１と、樹木頂点候補データが記憶されるファイル１２２と、樹冠上部抽出部１２３と、樹冠上部データが記憶されるファイル１２４と、樹木頂点決定部１２５と、樹木頂点データが記憶されるファイル１２６と出力部１３０等を備えている。
【００８４】
データベース１１０は、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model：本発明の第２の実施の形態では０．５ｍ間隔の格子状のデータ）を記憶している。このグリッドデータは、地盤のＸ、Ｙ、Ｚ（地盤高Ｚａ）とを対応させて記憶させられている。つまり、レーザーデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて地面を復元する。
【００８５】
データベース１１１は、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model：本発明の第２の実施の形態では０．５ｍ間隔）を記憶している。このグリッドデータ（本発明の第２の実施の形態では０．５ｍ間隔の格子状のデータ）は、樹木表面のＸ、Ｙ、Ｚ（Ｚｂ：樹木表面高）とを対応させて記憶している。つまり、レーザーデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて表層を復元する。
【００８６】
すなわち、これらのグリッドデータは、格子間隔０．５ｍのＸ−Ｙ平面直角座標系に、緯度経度に対応するＸ値、Ｙ値とレーザー計測によって得られたＺ値を割り当てている。
【００８７】
樹木高算出部１１４は、データベース１１０のＤＳＭとデータベース１１１のＤＥＭとの差を求め、これをファイル１１５（メモリ）のＸ，Ｙ座標に割付ける。
【００８８】
平滑フィルタ部１１６は、Average filterとも呼ばれるものであり、ファイル１１５の樹木高グリッドデータＤｉに対して、各々のデータを中心とするＸ−Ｙ座標面３×３範囲（本発明の第２の実施の形態では１．５ｍ×１．５ｍ）で平均を求め、これをファイル１２６に保存する。
【００８９】
尾根谷度算出部１１７は、ファイル１２６に平均化された樹木高グリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２・・）が保存されると、樹木高グリッドデータＤｄｉが付加されている任意の平均化された樹木高グリッドデータＤｄｉ（以下樹木高グリッドデータＤｄｉという）から各々のデータを中心としたＸ−Ｙ座標面の任意半径の検索範囲（検索範囲は樹冠の立体形状の特徴を適切に表現できた2mに決定している。以下２ｍ範囲毎）で尾根谷度を求め、これをファイル１１８に保存する。
【００９０】
局所最大値フィルタ部１２０は、ファイル１１８の尾根谷度のデータＤｅｉに３×３範囲のフィルタを順次適用し、この結果をファイル１２７に保存する。
【００９１】
樹木頂点候補抽出部１２１は、ファイル１２７の樹冠立体形状の特徴データＤｆｉ（ＸＹ、尾根谷度）とファイル１１８のＤｅｉ（Ｘｐｉ、Ｙｐｉ、尾根谷度）とを比較し、同じものを頂点候補Ｄｇｉとしてファイル１２２に保存する。
【００９２】
樹冠上部抽出部１２３は、ファイル１１８の尾根谷度データ（樹冠立体形状の特徴データＤｅｉ）に森林の状態に応じた閾値を設けることによって樹冠の上部のデータＤｅｉを抽出し、これを樹冠上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル１２４に保存する。
【００９３】
樹木頂点決定部１２５は、ファイル１２４の樹冠の上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル１２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定し、これをファイル１２８に保存する。そして、樹冠上部に複数の仮の頂点が存在するときは、この樹冠上部で最も樹木高が高いものを単木の頂点として決定し、これをファイル１２９に保存する（Ｐｉ）。
【００９４】
表示処理部１１９は決定した頂点を例えば○印（△印、×印、番号等の識別）にして画面に表示する。このとき、写真画像と重ねて表示してもよい。
【００９５】
出力部１３０は、ファイル１２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。このとき、樹木高も表示してもよい。
【００９６】
（動作説明）
上記のように構成された樹木頂点特定装置について以下に動作を説明する。本発明の第２の実施の形態では表示部１３にデータベース１１０のＤＳＭ及びデータベース１１１のＤＥＭ並びにデジタル写真を重ね表示している。
【００９７】
図１３は本発明の第２の実施の形態の樹木頂点認識装置の概略動作を説明するフローチャートである。
【００９８】
本発明の第２の実施の形態では、樹木高算出部１１４が航空レーザー計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ７１）。また、地盤のグリッドデータ（以下、ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ７２）。
【００９９】
そして、樹木高算出部１１４がDSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）をとり樹木高グリッドデータＤｉとして、ファイル１１５に保存する（Ｓ７３）。
【０１００】
この樹木高のグリッドデータＤｉに対して、平滑フィルタ部１１６（Average filter16ともいう）が３×３範囲で順次平均化し、ファイル１２６に保存（樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・））する（Ｓ７４）。
【０１０１】
次に、尾根谷度算出部１１７が２ｍ範囲毎（Ｄｄｉ）に地上開度及び地下開度を算出する（Ｓ７５）。次に、地下開度及び地上開度を用いて尾根谷度算出部１１７が「尾根谷度」データを作成することによって樹冠の立体形状の特徴を判別できるようにした尾根谷画像を得てファイル１１８に保存する（Ｓ７６）。
【０１０２】
尾根谷度は次の様に算出する。樹木高のグリッドデータＤｄｉを用いて、８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。
【０１０３】
検索範囲は1m・2m・5m・10mの4通り試算し、樹冠形状を適切に表現できた2mに決定した。
【０１０４】
そして、検索範囲2mにおける大まかな地形を表す、地上開度と地下開度に挟まれる角度の2等分線を尾根谷度（φ3）と定義していく（図２０参照）。この尾根谷度のデータＤｅｉは表示処理部１１９によって表示（尾根谷度画像）される。
【０１０５】
さらに、尾根谷度のデータＤｅｉに局所最大値フィルタ部１２０を順次適用し（Ｓ７７）、これをファイル１２７に保存する。
【０１０６】
このとき、樹冠の大きさに限らず最小のエリアサイズである3×3範囲に固定して、樹頂点ではなく「樹頂点の候補点」を抽出し、ファイル１２２に樹木頂点候補Ｄｇｈ（ｄｇ１、・・）として保存する（Ｓ７８）。この３×３範囲に固定したのは、例えば５×５にすると、樹木頂点候補が存在しない樹冠上部が多くなるという実験結果を得たので、３×３範囲とした。
【０１０７】
一方で、「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった性質を利用して、樹木頂点候補抽出部１２１が尾根谷度データに森林の状態に応じた閾値を設けることによって樹冠の上部を抽出して（Ｓ７９）、冠上部抽出部１２４をラベリング化する（Ｓ８０）。
【０１０８】
さらに、樹木頂点決定部１２５が樹頂候補点と樹冠上部エリア毎に、樹木頂点候補点の中から最高樹木高を求め、その点を樹木頂点と決定（Ｓ８１）し、これを樹木頂点として抽出する（Ｓ８２）。この樹木頂点の決定については後述する。
【０１０９】
図１４は本発明の第２の実施の形態の樹木頂点認識装置の尾根谷度を求める処理を説明するフローチャートである。
【０１１０】
本発明の第２の実施の形態では、表示部１３にオルソフォト画像を表示している（Ｓ９０）。
【０１１１】
樹木高算出部１１４は、ＤＳＭ（Ｘ、Ｙ、Ｚ）及びＤＥＭ（Ｘ、Ｙ、Ｚ）とを読み込み（Ｓ９１）、差Ｚｃｉを求め（Ｓ９２）、これらを樹木高のグリッドデータＤｉ（Ｄ１、Ｄ２、・・・）としてファイル１１５に保存する。
【０１１２】
この樹木高データＤｉは、図１５に示すように、座標（Ｘ，Ｙ：平面直角座標）と樹木高Ｚｃｉとが対応させられたデータとなる。
【０１１３】
次に、平滑フィルタ部１１６がファイル１１５の多数の樹木高のグリッドデータＤｉに対してフィルタ１６で順次平均化（３×３範囲）し（Ｓ９３）、この樹木高のグリッドデータＤｄｉをファイル１２６に保存する。
【０１１４】
次に、尾根谷度算出部１１７は、個々のDdiにおいて２ｍ範囲で順次、尾根谷度を求めることによって樹冠の立体形状の特徴を判別する（Ｓ９４）。
【０１１５】
尾根谷度は、８方向に探索（木の種類によって１ｍ、５ｍ、１０ｍで行う場合もある）して地上開度と地下開度とを求め、地上開度と地下開度との差を２で割ることで求める（Ｓ９４）。つまり、樹木の冠部（葉、枝）の立体形状の特徴が求められる。
【０１１６】
ファイル１２６の樹木高のグリッドデータＤｄｉを用いて、図１６に示すように８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。検索範囲は樹冠の立体形状の特徴を適切に表現できた2mに決定している。
【０１１７】
図１７は任意の検索範囲における樹木高のグリッドデータＤｄｉを距離と高さからなる軸にプロットしたものである。また、図１７の黒点は０．５ｍ間隔のグリッドデータＤｄｉを示す（２ｍ範囲は図示せず）。
【０１１８】
つまり、図１７は着目する標本地点から任意半径の範囲で一方向を見た場合の大まかな樹冠形状を表している。そして、この図１７において、φ１が地上開度（地上角）であり、φ２が地下開度（地下角）である。そして、図１８（ａ）（ｂ）では地上開度と地下開度に挟まれる角度の2等分線と距離軸（水平軸）とで挟まれた角度を尾根谷度（φ3）と定義している。また、図１８の黒点は０．５ｍ間隔のグリッドデータＤｄｉを示す（２ｍ範囲は図示せず）
本形態では索範囲を２ｍとして、尾根谷度（φ3）を求め、尾根谷度のデータ（以下樹冠立体形状の特徴データＤｅｉ（Ｄｅｉ：Ｄｅ１、Ｄｅ２・・・）としてファイル１１８に保存する（Ｓ９６）。
【０１１９】
図１９が樹木高をY軸とし距離をＸ軸として軌跡Ａｉの樹冠形状を表現した説明図である。図２０が尾根谷度をY軸にし距離をＸ軸として軌跡Ａｉの樹冠の立体形状の特徴を表現した説明図である。
【０１２０】
つまり、図２０に示すように尾根谷度で軌跡Ａｉの樹冠の立体形状の特徴を表現すると、図１９より顕著に樹冠間の境界が現れる。図１９及び図２０の番号（１）〜（６）は写真から肉眼で読み取った樹冠エリアを示すものである。
【０１２１】
単木としての認識しやすさを向上するには、樹頂点とエッジ部分の値の差が大きいこと、樹頂点において値が同傾向を示すこと、エッジにおいて値が同傾向を示すこと等が挙げられる。図２０に示すように、従来までの「高さ」データより本手法の「尾根谷度」データの方が先にあげた条件を満たしていることがわかる。これにより、樹冠ポリゴンを自動作成する際の精度向上に寄与することができる。
【０１２２】
図２０においては、尾根谷度がマイナスの値となる点から急激に＋方向にあがり、そして急激にマイナスの値となる点までの範囲が単木の樹冠部であることが分かる。
【０１２３】
また、図２２は頂点候補と樹頂点の選定法を説明する図であり、各樹冠上部で樹木高が最も高い候補を選定することし示している。この樹木頂点の決定については詳細に後述する。
【０１２４】
次に、尾根谷度算出部１１７は、求めた尾根谷度を樹冠立体形状の特徴データＤｅｉとしてファイル１１８（ＸＹが対応させられている）に保存（図２１参照）する（Ｓ９６）。
【０１２５】
前述の地上開度は、各要素（Ｄｄｉ）に対し、その要素から8方向に探索し、考慮距離内で最大傾斜角を求める。8方向それぞれの最大傾斜角の平均を、天頂からの角度で表す。なお、地上開度の値は考慮範囲に依存し、判読を行いたい対象のスケールに合わせて設定する必要がある。
【０１２６】
また、地下開度は各要素に対し、その要素から8方向に探索し、考慮距離内で最小傾斜角を求める。8方向それぞれの最小傾斜角の平均を、天頂からの角度で表す。地上開度と同様に、値は考慮範囲に依存し、判読を行いたい対象のスケールに合わせて設定する必要がある。
【０１２７】
尾根谷度は、算出された地上開度から地下開度を引き、２で割ることによって計算する。尾根地形は値がプラスになり、谷地形は値がマイナスで表現される。
【０１２８】
そして、図１３の樹木頂点抽出処理に進む。樹木頂点抽出処理のステップＳ７でファイル１１８の樹冠立体形状の特徴データＤｅｉ（Ｄｅ１、Ｄｅ２・・）に対して３×３の局所最大値フィルタ部１２０（３×３範囲のウインドウを用いる）を適用して、これをファイル１２７に保存（Ｄｆｉ）する。
【０１２９】
ステップＳ７８で、樹木頂点候補抽出部１２１が樹冠エリア（(1)、(2)・・）に存在するピーク点の樹冠立体形状の特徴データＤｆｉ（図２２参照）を頂点候補（以下Ｄｇｉ）として抽出し、ファイル１２２に保存する（図２３参照）。具体的にはファイル１１８のデータの値とファイル１２７の同じ座標値のデータの値が同値の場合を「１」と、他を「０」としてファイル１２２に定義する（図２３参照）。
【０１３０】
ここで、局所最大フィルタ部１２０及び樹木頂点候補抽出部１２１の処理について説明を補充する。
【０１３１】
つまり、局所最大フィルタ部１２０（Local max filter２０ともいう）は、ファイル１１８の樹冠立体形状の特徴データＤｅｉ（Ｄｅ１、Ｄｅ２・・）に３×３範囲のウインドウをかけ、このウインドウ内で最大となる値をウィンドウの中央値に置き換える処理を順次行う。これらを樹冠立体形状の特徴データＤｆｉとしてファイル１２７に保存する。
【０１３２】
そして、樹木頂点候補算出部２１が樹冠立体形状の特徴データＤｆｉ（ＸＹ、尾根谷度）とファイル１１８のＤｅｉ（ＸＹ、尾根谷度）とを比較し、同じものを頂点候補Ｄｇｉとしてファイル１２２に保存している。
【０１３３】
一方、ステップＳ７９で樹冠上部抽出部１２３が「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった性質を利用して、ファイル１１８の尾根谷度データ（樹冠立体形状の特徴データＤｅｉ）に森林の状態に応じた閾値を設けることによって樹冠の上部のデータＤｅｉを抽出し、これを樹冠上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル１２４に保存する（図２４参照）。
【０１３４】
図２４（ａ）の樹冠の上部のデータＤｈｉは、閾値を超えるＤｅｉを抽出したものであり、例えば番号１のエリアのデータＤｅｉは、（Ｘｈ５、Ｙｈ５）を「１」、（Ｘｈ６、Ｙｈ６）を「１」・・（Ｘｈ１０、Ｙｈ１０）を「１」として示している。
【０１３５】
図２４の（ｂ）は図２４（ａ）を具体的に示すものであり、閾値を超えた頂点は「１」とする。「１」が連続している範囲が樹冠上部である。
【０１３６】
そして、ステップＳ８０でファイル１２４の中の頂点候補ＤｈｉのＸｈｉ、Ｙｈｉをポリゴン化（座標Ｘｈｉ、Ｙｈｉと色情報）して表示処理部１１９によってラベルリング（纏まったエリアを（上部）にラベルを振る）させる（図２５参照）。
【０１３７】
次に、ステップＳ８１で樹木頂点決定部１２５は、ファイル１２４の樹冠の上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル１２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定する。そして、樹冠上部のエリア内に複数の仮の頂点が存在するときは、この樹冠上部のエリア（図１９、図２０においては、番号１、２・・のエリア）で最も樹木高が高いものを単木の頂点として決定し、これをファイル１２９に保存する（Ｐｉ）。表示処理部１１９は決定した頂点を画面に表示する（図２６参照）。すなわち、例えば、頂点候補Ｄｇｉは図２７に示すように表示されるが、頂点決定された画像は図２８に示すようになる。
【０１３８】
前述の頂点の決定について説明を補充する。例えば、図２９の（ａ）に示すように、ファイル１２４に第１列目が「００１１００１１１００」、第２列目が「０１１１１０１１１１０」、第３列目が「００１１０００１０００」、第４列目が「０００００１０００１１」という具合に頂点上部域が得られた場合に、図２９（ｂ）に示すように、「１」が隣接するまとまり（「０」が区切り）を１つの樹冠上部としてまとまり毎に通し番号を割り振っていく。
【０１３９】
そして、例えば、同一高の頂点が樹冠内に複数の頂点が存在することになると、一度、ファイル１２８のデータＤＱｉを表示させて、オペレータによって指定させ、指定されたものを頂点Ｐｉと決定してファイル１２９に保存する。
【０１４０】
次に、尾根谷度について説明を補充する。尾根谷度の特徴としては、図３０に示すようなsinカーブの地形を考えると、尾根地形は「尾根谷度＞0」、谷地形は「尾根谷度＜0」で表現されるとともに、中腹のところで「尾根谷度＝0」となることが挙げられる。
【０１４１】
ここで、わかりやすくするために、理想の樹木が図３０に示すような形状をしているとすると、樹頂点では「尾根谷度＞0」となっており、必ず「尾根谷度＝0」となる点を通過して樹冠のエッジに到達し、そこでは「尾根谷度＜0」となっていると考えられる。したがって、樹冠の上部は「尾根谷度＞0」で定義できると言える。
【０１４２】
このことを利用して、本業務では、図３１に示すように樹頂点を絞り込むこととした。つまり、尾根谷度の値を利用することによって樹冠上部を抽出し、この抽出した中に存在する樹頂候補点のうち、最も樹高（DSM-DEM）が高い候補点を樹頂点として選定することとした。
【０１４３】
尾根谷度の閾値について説明を補充する。
【０１４４】
尾根谷度に閾値を設定して樹冠上部を抽出する際、図３２に示すように閾値αの場合には樹木５が抽出できなくなってしまうことがわかる。一方で、閾値βの場合には樹木５と樹木６が一本として認識されてしまう。このように、閾値によって抽出される樹頂点が左右されしまう問題を抱えている。
【０１４５】
したがって、ここでは、樹冠上部を抽出するための最適な閾値について検討することとした。具体的には、尾根谷度の閾値を「10,15,20,25,30」の５種類に設定して、樹頂点の抽出結果が最も良いものを採用することとした。絞込み後の樹頂点の集計結果を表１に示す。
【表１】
【０１４６】
表１を見ると、「尾根谷度＞20」での結果が最もよかったことがわかる。これと比較して「尾根谷度＞10」での抽出本数が少ない理由は隣接する樹冠が一つとして認識されたためである。また、「尾根谷度＞30」での抽出本数が少ない理由は、閾値を高く設定したために抽出できなかった樹冠があったためである。
【０１４７】
従って、樹冠上部を抽出するための尾根谷度の閾値は「20」が最適であることがわかった。樹冠上部の抽出画像は図２５、絞込み後の樹頂点抽出画像は図２６に示している。
【０１４８】
上記のようにある山林の樹冠の頂点を算出することによって、単木を特定できるので山林の資産の管理等に適用できる。例えば、出力部１３０がファイル１２９の頂点のデータＰｉ（Ｘ、Ｙ、Ｚ）を出力して印刷させると、本数と高さが分かるので山林の価値等が自動的に把握できる。
【０１４９】
図３３においては、樹木頂点を算出した後で樹冠ポリゴンを作成している。樹冠ポリゴン作成は、ステップＳ８２で抽出された樹木頂点とファイル１１８の尾根谷度を用いてWetershedアリゴリズムにより、樹冠ポリゴン（図３４参照）を生成する（Ｓ８３、Ｓ８５）。
【０１５０】
また、ステップＳ８３の閾値の決定方法は下記の通りである。
【０１５１】
図１４の写真画像のようにオルソフォト画像において任意の線（軌跡Ａｉ）を引き、尾根谷度算出部が尾根谷度の横断図を作成し、これを表示する（図２０に該当）。
【０１５２】
次に、この横断線のスタートからの追加距離（図２０のX軸に相当）を追い、表示部のオルソ画像の樹冠と樹冠の境界を探す。
【０１５３】
そして、境界位置の追加距離を見て、その追加距離の尾根谷度（尾根谷度算出部が算出して表示させる）を図２０から作業者が読み取る。
【０１５４】
この行程を繰り返し、樹冠上部が形成される尾根谷度を読み取る。
【０１５５】
図２０の(1)では0と-50、(2)では-50と-20、(3)では-20と0、(4)では0と0、(5)では0と10、(6)では10と-10となった。
【０１５６】
この結果から、(1)〜(6)が分かれる尾根谷度10を閾値とする。
【０１５７】
閾値を0とすると(5)と(6)が一つの樹冠上部になるため、0は不採用とする。
【０１５８】
閾値を30とすると(2)では樹冠上部が2つ形成され、(5)が樹冠上部にならないので、これも不採用とする。よって、図２０では尾根谷度10程度が閾値としてふさわしいことになる。このようにして決定した閾値を用いている。
【０１５９】
なお、上記実施の形態では航空写真を用いて説明したが、樹木高のグリッドデータ又はＤＥＭ、ＤＳＭにグレイスケールを割り当ててもよい。
【０１６０】
例えば、地上開度と地下開度の差分画像をグレイに、傾斜を赤のチャンネルにいれて、擬似カラー画像を作成することにより、尾根や山頂部分を白っぽく、また谷や窪地を黒っぽく表現し、傾斜が急な部分ほど赤く表現する。このような表現の組み合わせにより、１枚でも立体感のある画像（立体赤色化マップともいう）を用いてもよい。
【０１６１】
例えば、地上開度の値が４０度から１２０度程度の範囲に収まる場合は、５０度から１１０度をグレイスケールに対応させ、２５５諧調に割り当てる。つまり、尾根の部分（凸部）の部分ほど地上開度の値が大きいので、色が白くなる。
【０１６２】
また、本発明の第２の実施の形態は、図３５に示すようなコンピュータシステムVPSによって実現するのが好ましい。ＶＰＳ１は、ワークステーション、プロセッサ、マイクロコンピュータ、ロジック、レジスタ等の適宜な組み合わせからなる中央情報処理装置（CPU）と、この中央情報処理装置に必要な制御・操作情報を入力するキーボード（KB）、マウス、対話型ソフトスイッチ、外部通信チャネル等を含む情報入力部と、中央情報処理装置から出力された情報を広義な意味で表示・伝送するディスプレイ、プリンタ、外部通信チャネル等を含む情報出力部と、中央情報処理装置に読み込まれるオペレーティングシステム、アプリケーションプログラム等の情報が格納されたロム（ＲＯＭ）と、中央情報処理装置で随時処理すべき情報及び中央情報処理装置から随時書き込まれる情報を格納するラム（ＲＡＭ）等とを備える。ＲＯＭ、ＲＡＭを適宜統合、細分化することは差し支えない。
【０１６３】
前述のアプリケーションプログラムは、本発明の第２の実施の形態の、樹木高算出部１１４と、尾根谷度算出部１１７と、表示処理部１１９と、局所最大フィルタ部１２０（Ｌｏｃａｌ max ｆｉｌｔｅｒともいう）と、樹木頂点候補抽出部１２１と、樹冠上部抽出部１２３と、樹木頂点決定部１２５と、出力部１３０等である。
【０１６４】
本実施形態では、開度という概念を用いている。開度は当該地点が周囲に比べて地上に突き出ている程度及び地下に食い込んでいる程度を数量化したものである。つまり、地上開度は、着目する標本地点から距離Ｌの範囲内で見える空の広さを表しており、また地下開度は逆立ちをして地中を見渡す時、距離Ｌの範囲における地下の広さを表している。
【０１６５】
開度は距離Ｌと周辺地形に依存している。一般に地上開度は周囲から高く突き出ている地点ほど大きくなり、山頂や尾根では大きな値をとり窪地や谷底では小さい。逆に地下開度は地下に低く食い込んでいる地点ほど大きくなり、窪地や谷底では大きな値をとり山頂や尾根では小さい。
【０１６６】
また、開度図は計算距離の指定によって、地形規模に適合した情報抽出が可能であり、方向性及び局所ノイズに依存しない表示が可能である。
【０１６７】
つまり、尾根線及び谷線の抽出に優れており、豊富な地形・地質情報が判読できるものであり、一定範囲のＤＥＭデータ上において、設定した当該地点Ａから８方向のいずれか一方を見たときに最大頂点となる点Ｂを結ぶ直線Ｌ１と、水平線とがなす角度ベクトルθｉを求める。この角度ベクトルの求め方を８方向に渡って実施し、これらを平均化したものを地上開度と称し、一定範囲のＤＥＭデータ上（地表面：立体）に空気層を押し当てた立体を裏返した反転ＤＥＭデータの当該地点Ａから８方向のいずれか一方を見たときに最大頂点となる点Ｃ（一番深い所に相当する）を結ぶ直線Ｌ２と、水平線とがなす角度を求める。この角度を８方向に渡って求めて平均化したのを地下開度ψｉと称している。
【０１６８】
すなわち、地上開度は、着目点から一定距離までの範囲に含まれるＤＥＭデータ上において、８方向毎に地形断面を生成し、それぞれの地点と着目点を結ぶ線Ｌ１の傾斜の最大値（鉛直方向から見たとき）を求める。このような処理を８方向に対して行う。傾斜の角度は天頂からの角度（平坦なら９０度、尾根や山頂では９０度以上、谷底や窪地では９０度以下） また、地下開度は、反転ＤＥＭデータの着目点から一定距離までの範囲において、８方向毎に地形断面を生成し、それぞれの地点と着目点を結ぶ線の傾斜の最大値（地表面の立体図において鉛直方向からＬ２を見たときには最小値）を求める。このような処理を８方向に対して行う。
【０１６９】
地表面の立体図において鉛直方向からＬ２を見たときの角度ψｉは、平坦なら９０度、尾根や山頂では９０度以下、谷底や窪地では９０度以上である。
【０１７０】
つまり、地上開度と地下開度は、２つの基本地点Ａ（ｉA，ｊA，ＨA）とＢ（ｉB，ｊB，ＨB）を考える。
【０１７１】
ＡとＢの距離は
P = {(iA − iB)2 + (jA − jB)2}1/2 …(1)
となる。
【０１７２】
標高０ｍを基準として、標本地点のＡとＢの関係を示したものである。標本地点Ａから標本地点Ｂを見た仰角θは
θ＝ｔａｎ-1｛（ＨB−ＨA）／Ｐ｝
で与えられる。θの符号は(1)ＨA＜ＨB の場合には正となり、(2)ＨA＞ＨB の場合には負となる。
【０１７３】
着目する標本地点から方位Ｄ距離Ｌの範囲内にある標本地点の集合を DＳL と記述して、これを「着目する標本地点のＤ−Ｌ集合」を呼ぶことにする。ここで、
DβL ：着目する標本地点の DＳL の各要素に対する仰角のうちの最大値
DδL ：着目する標本地点の DＳL の各要素に対する仰角のうちの最小値
として、次の定義をおこなう。
【０１７４】
定義１：着目する標本地点のＤ−Ｌ集合の地上角及び地下角とは、各々
DφL ＝９０− DβL
及び
DψL ＝９０＋ DδL
を意味するものとする。
【０１７５】
DφL は着目する標本地点から距離Ｌ以内で方位Ｄの空を見ることができる天頂角の最大値を意味している。一般に言われる地平線角とはＬを無限大にした場合の地上角に相当している。また、DψL は着目する標本地点から距離Ｌ以内で方位Ｄの地中を見ることができる天底角の最大値を意味している。Ｌを増大させると、 DＳL に属する標本地点の数は増加することから、 DβL に対して非減少特性を持ち、逆に DδL は非増加特性を持つ。したがって DφL 及びDψLは共にＬに対して非増加特性を持つことになる。
【０１７６】
測量学における高角度とは、着目する標本地点を通過する水平面を基準にして定義される概念であり、θとは厳密には一致しない。また地上角及び地下角を厳密に議論しようとすれば、地球の曲率も考慮しなければならず、定義１は必ずしも正確な記述ではない。定義１はあくまでもＤＥＭを用いて地形解析をおこなうことを前提として定義された概念である。
【０１７７】
地上角及び地下角は指定された方位Ｄについての概念であったが、これを拡張したものとして、次の定義を導入する。
【０１７８】
定義II：着目する標本地点の距離Ｌの地上開度及び地下開度とは、各々 ΦL＝（0φL ＋45φL ＋90φL ＋135φL ＋180φL ＋225φL ＋270φL ＋315φL ）／８
及び
ΨL＝（０ψL ＋45ψL ＋90ψL ＋135ψL ＋180ψL ＋225ψL ＋270ψL ＋315ψL ）／８
を意味するものとする。
【０１７９】
地上開度は着目する標本地点から距離Ｌの範囲内で見える空の広さを表しており、また地下開度は逆立ちをして地中を見渡す時、距離Ｌの範囲における地下の広さを表している。
【０１８０】
（その他の実施の形態）
上記のように、本発明は第１及び第２の実施の形態によって記載したが、この開示の一部をなす論述及び図面はこの発明を限定するものであると理解すべきではない。この開示から当業者には様々な代替実施の形態、実施例及び運用技術が明らかとなろう。
【０１８１】
例えば、本発明の第１の実施の形態では、樹冠反射強度抽出部１３が、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ifirstを、樹冠反射強度Ictとして抽出する例を説明したが、樹木頂点Ｐiの代わりに、本発明の第２の実施の形態において説明した樹冠上部データ抽出部１２３により抽出された樹冠上部のデータＤhiを用いても良い。即ち、樹冠反射強度抽出部１３が、樹冠上部のデータＤhiを用いて、樹冠上部付近の座標（Ｘ，Ｙ）のラストパルスの強度Ifirstを、樹冠反射強度Ictとして抽出しても良い。
【０１８２】
また、本発明の第１の実施の形態では、樹冠通過強度抽出部１４が、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ilastを、樹冠通過強度Icbとして抽出する例を説明したが、樹木頂点Ｐiの代わりに、本発明の第２の実施の形態において説明した樹冠上部データ抽出部１２３により抽出された樹冠上部のデータＤhiを用いても良い。即ち、冠通過強度抽出部１４が、樹冠上部のデータＤhiを用いて、樹冠上部付近の座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ilastを、樹冠通過強度Icbとして抽出しても良い。
【０１８３】
また、本発明の第２の実施の形態において樹木頂点の決定方法を説明したが、本発明の第１の実施の形態において用いる樹木頂点は、本発明の第２の実施の形態において説明した決定方法によるものでなくても良く、決定方法は特に限定されない。例えば、樹木頂点の特定方法としては、写真画像等の樹木の頂点を目視により特定する手法や、航空レーザー計測により取得されたデータを用いてＤＳＭとＤＥＭを作成し、ＤＳＭとＤＥＭの差分をとることによって樹木高のグリッドデータを算出し、樹木高や樹木高の極大値を判別指標として樹木頂点を特定する手法を用いても良い。
【０１８４】
このように、本発明はここでは記載していない様々な実施の形態等を含むことは勿論である。したがって、本発明の技術的範囲は上記の説明から妥当な特許請求の範囲に係る発明特定事項によってのみ定められるものである。
【図面の簡単な説明】
【０１８５】
【図１】本発明の第１の実施の形態に係る具体的なハードウエア構成図である。
【図２】本発明の第１の実施の形態に係る葉面積指数算出装置の概略構成図である。
【図３】レーザー反射強度を説明する説明図である。
【図４】分光特性（反射率）のグラフである。
【図５】算出エリアを説明するための平面図である。
【図６】抽出した樹冠内部情報（レーザー反射強度）を示す表である。
【図７】算出した葉面積指数を示す表である。
【図８】本発明の第１の実施の形態に係る葉面積指数算出方法のフローチャートである。
【図９】地盤反射強度の設定方法のフローチャートである。
【図１０】樹冠反射強度の設定方法のフローチャートである。
【図１１】樹冠通過強度の設定方法のフローチャートである。
【図１２】本発明の第２の実施の形態の樹木頂点認識装置の概略構成図である。
【図１３】本発明の第２の実施の形態の樹木頂点認識装置の概略動作を説明するフローチャートである。
【図１４】本発明の第２の実施の形態の樹木頂点認識装置の尾根谷度を求める処理を説明するフローチャートである。
【図１５】樹木高データＤiの説明図である。
【図１６】尾根谷度の算出を説明する説明図である。
【図１７】尾根谷度の算出を説明する説明図である。
【図１８】尾根谷度の算出を説明する説明図である。
【図１９】樹木高をＹ軸とし距離をＸ軸として軌跡Ａｉの樹冠形状を表現した説明図である。
【図２０】尾根谷度をＹ軸にし距離をＸ軸として軌跡Ａｉの樹冠の立体形状の特徴を表現した説明図である。
【図２１】樹木高グリッドデータの説明図である。
【図２２】尾根谷度による頂点候補の算出を説明する説明図である。
【図２３】頂点候補Ｄgiの説明図である。
【図２４】樹冠の上部のデータＤhiの説明図である。
【図２５】樹冠上部のラベリングを説明する説明図である。
【図２６】決定した頂点を画面の説明図である。
【図２７】頂点決定までの説明図である。
【図２８】頂点決定までの説明図である。
【図２９】頂点決定までの説明図である。
【図３０】尾根谷度の説明図である。
【図３１】尾根谷度の説明図である。
【図３２】閾値と樹冠部の関係を説明する説明図である。
【図３３】他の実施の形態を説明するフローチャートである。
【図３４】他の実施の形態によって得られた樹冠のポリゴンを説明する説明図である。
【図３５】本発明の第２の実施の形態の具体的なハードウエア構成図である。
【符号の説明】
【０１８６】
１…ＣＰＵ
２…表示部
３…入力部
４…出力部
５…記憶装置
１０…表示処理部
１１…地盤反射強度抽出部
１２…I0算出部
１３…樹冠反射強度抽出部
１４…樹冠通過強度抽出部
１５…葉面積指数算出部
１６…平滑フィルタ（アベレージフィルタ）
２０…オルソフォト画像記憶部
２１…反射強度データ記憶部
２２…反射率記憶部
２３…樹木頂点記憶部
２４…ＤＥＭ記憶部
２５…閾値記憶部
２６…地盤反射強度記憶部
２７…入射強度記憶部
２８…樹冠反射強度記憶部
２９…樹冠通過強度記憶部
３１…葉面積指数記憶部
１１０…データベース
１１１…データベース
１１２…データベース
１１３…表示部
１１４…樹木高算出部
１１６…平滑フィルタ
１１７…尾根谷度算出部
１１８…ファイル
１１９…表示処理部
１２０…局所最大値フィルタ
１２１…樹木頂点候補抽出部
１２２…ファイル
１２３…樹冠上部データ抽出部
１２４…ファイル
１２５…樹木頂点決定部
１２６…ファイル
１２７…ファイル
１２８…ファイル
１２９…ファイル
１３０…出力部
【技術分野】
【０００１】
本発明は、林分パラメータとして重要な葉面積指数を算出する葉面積指数算出装置、葉面積指数算出方法及びそのプログラムに関する。
【背景技術】
【０００２】
樹木本数、樹冠面積、樹高、胸高直径、及び葉面積指数(Leaf Area Index:ＬＡＩ)等の林分パラメータを把握することは、森林の状態を定量化し、林業や保全活動へ反映させる上で重要である。ここで、ＬＡＩは、植物の葉の多少の度合いを示す指標であり、地上の単位面積に対しての、その上方に存在するすべての葉の片側の総面積の比率である（例えば、特許文献１参照。）。
【０００３】
ＬＡＩ推定手法としては、現地調査を行うのが一般的である。しかし、この方法では、局所的にしか推定することができず、広域のＬＡＩを推定するのには、時間や労力がかかる。
【０００４】
そこで、広域にＬＡＩを推定する場合には、衛星データを活用する手法がある。これは、衛星データより正規化植生指標（ＮＤＶＩ）を求め、ＮＤＶＩとＬＡＩの回帰分析を実施することによって得られた回帰式に基づいて推定をするものである。この手法においても、ＬＡＩの実測値を取得していることが前提であり、広域展開できるものの、実測データがなければ実施できない。
【０００５】
また、航空レーザーを活用して広域のＬＡＩを推定する手法がある。この方法では、航空レーザーにより計測されたレーザー反射地点の高さとＤＥＭの差をとった地上面からの高さ（樹高データ）を用いてＬＡＩを推定する。具体的には、横軸を地上面からの高さ、縦軸を頻度としたヒストグラムから累積頻度分布を作成し、これが光の減衰を表現していると仮定してＬＡＩを算出している。
【特許文献１】特開２００５−５２０４５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来の航空レーザーを活用したＬＡＩ推定手法では、ヒストグラムが樹冠内部の情報を反映しきれていないため、算出したＬＡＩには疑問が残る。即ち、ＬＡＩを高精度に推定することができない問題点があった。
【０００７】
上記問題点を鑑み、本発明は、航空レーザー計測を用いて広域かつ簡易にＬＡＩを高精度に推定することができる葉面積指数算出装置、葉面積指数算出方法及びそのプログラムを提供することである。
【課題を解決するための手段】
【０００８】
本願発明の一態様によれば、（イ）上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶するステップと、（ロ）地盤の反射率を記憶するステップと、（ハ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出するステップと、（ニ）記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出するステップと、（ホ）記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出するステップと、（ヘ）記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出するステップと、（ト）樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出するステップとを含む葉面積指数算出方法が提供される。
【０００９】
本願発明の他の態様によれば、コンピュータを、（イ）上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手段と、（ロ）地盤の反射率を記憶する手段と、（ハ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手段と、（ニ）記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手段と、（ホ）記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出する手段と、（ヘ）記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出する手順と、（ト）樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する手段手段として機能させるためのプログラムが提供される。
【００１０】
本願発明の更に他の態様によれば、（イ）上空から地上へレーザーを照射するレーザー照射手段と、（ロ）地上で反射したレーザーの反射強度データを記憶する反射強度記憶手段と、（ハ）地盤の反射率を記憶する反射率記憶手段と、（ニ）記憶した反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出し、記憶した反射強度データから、樹冠を通過して地盤で反射したレーザーの反射強度を抽出し、記憶した反射強度データから、地盤で反射したレーザーの反射強度を抽出し、記憶した反射強度データ及び反射率に基づいて、樹冠に到達する直前のレーザーの入射強度を算出し、樹冠に到達する直前のレーザーの入射強度、樹冠を通過したレーザーの反射強度、樹冠で反射したレーザーの反射強度、及び地盤で反射したレーザーの反射強度を用いて、樹木の葉面積指数を算出する演算手段と、（ホ）葉面積指数を出力する出力手段とを備える葉面積算出装置が提供される。
【発明の効果】
【００１１】
本発明によれば、航空レーザー計測を用いてＬＡＩを高精度に推定することができる葉面積指数算出装置、葉面積指数算出方法及びそのプログラムを提供することができる。
【００１２】
したがって、このＬＡＩを林分パラメータとして用いて、ＣＯ2削減量を精度良く推定すること等が可能となる。
【発明を実施するための最良の形態】
【００１３】
次に、図面を参照して、本発明の第１及び第２の実施の形態を説明する。以下の図面の記載において、同一又は類似の部分には同一又は類似の符号を付している。また、以下に示す実施の形態は、この発明の技術的思想を具体化するための装置や方法を例示するものであって、この発明の技術的思想は、構成部品の材質、形状、構造、配置等を下記のものに特定するものでない。この発明の技術的思想は、特許請求の範囲において、種々の変更を加えることができる。
【００１４】
本発明の第１の実施の形態に係る葉面積指数算出装置は、図１に示すように、中央処理装置（ＣＰＵ）１、表示部２、入力部３、出力部４及び記憶装置５を備える。入力部３としては、例えばキーボード、マウス、ＯＣＲ等の認識装置、スキャナ、カメラ等の画像入力装置、マイク等の音声入力装置等が使用可能である。出力部４としては、液晶ディスプレイ、ＣＲＴディスプレイ等の表示装置や、インクジェットプリンタ、レーザープリンタ等の印刷装置等を用いることができる。表示部２は、液晶ディスプレイ、ＣＲＴディスプレイ等である。記憶装置５には、ＲＯＭ及びＲＡＭが組み込まれている。ＲＯＭは、ＣＰＵ１において実行されるプログラムを格納する。ＲＡＭは、ＣＰＵ１におけるプログラム実行処理中に利用されるデータ等を一時的に格納したり、作業領域として利用される一時的なデータメモリ等として機能する。記憶装置５としては、例えば半導体メモリ、磁気ディスク、光ディスク、光磁気ディスクや磁気テープ等が採用可能である。
【００１５】
図２は、本発明の第１の実施の形態の葉面積指数算出装置の概略構成図である。（第１の）地盤反射強度抽出部６、入射強度算出部７、表示処理部１０、（第２の）地盤反射強度抽出部１１、反射率算出部１２、樹冠反射強度抽出部１３、樹冠通過強度抽出部１４、及び葉面積指数算出部１５は、図１に示したＣＰＵ１により実現される。
【００１６】
図２に示す地盤反射強度記憶部８、反射率記憶部９、オルソフォト画像記憶部２０、反射強度データ記憶部２１、反射率記憶部２２、樹木頂点記憶部２３、ＤＥＭ記憶部２４、閾値記憶部２５、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８、樹冠通過強度記憶部２９及び葉面積指数記憶部３１は、図１に示した記憶装置５に含まれる。
【００１７】
オルソフォト画像記憶部２０は、ＬＡＩ算出対象となる森林域Ａを含むオルソフォト画像を格納する。「オルソフォト画像」とは、写真画像に三次元計測データ等を与えて正射変換を行った画像である。オルソフォト画像は、例えば、座標（Ｘ，Ｙ）とそれに対応する赤（Ｒ）、緑（Ｇ）、青（Ｂ）の色情報を含む。
【００１８】
反射強度データ記憶部２１は、航空レーザー（以下、単に「レーザー」という。）の反射強度データを格納する。反射強度データは、航空機に設けたレーザー照射手段（図示省略）によりレーザーを地上へ照射し、航空機内のレーザー受信手段（図示省略）により地上で反射したレーザーを受信することで取得することができる。反射強度データは、任意の地点で得られるランダムデータであり、３次元座標（Ｘ，Ｙ，Ｚｄ）とレーザーの反射強度を含む。高さの座標Ｚｄは、地盤や樹木等のレーザーの反射地点の高さを示す。
【００１９】
レーザーの反射強度は、パルス波形で示される。パルス波形には複数のピーク値があり、最初のピーク値を「ファーストパルスの強度Ifirst」、最後のピーク値を「ラストパルスの強度Ilast」と定義する。
【００２０】
ＤＥＭ記憶部２４は、森林の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を格納する。ＤＥＭは、例えば平面直角座標系のＸ軸、Ｙ軸を一定間隔（０．５ｍ）で区切ったメッシュの格子点に割り当てている。ＤＥＭは、３次元座標（Ｘ，Ｙ，Ｚh）を含む。高さの座標Ｚｈは、地盤の高さを示す。
【００２１】
地盤反射強度抽出部６は、反射強度データ記憶部２１から、既知反射率領域内の座標（Ｘ，Ｙ）の複数の反射強度データに含まれるファーストパルスの強度Ifirstをそれぞれ読み出す。ここで、「既知反射率領域」とは、地盤の反射率が予め分かっている領域をいい、例えばアスファルトや草原の領域である。更に、地盤反射強度抽出部６は、複数のファーストパルスの強度Ifirstの平均値を算出し、その平均値を、図３に示すような、既知反射率領域における地盤に直接達して反射したレーザーの反射強度（以下、「地盤反射強度」という。）Iaとして抽出する。既知反射率領域における地盤反射強度Iaは、地盤反射強度記憶部８に記憶される。
【００２２】
入射強度算出部７は、地盤反射強度記憶部８から、既知反射率領域における地盤反射強度Iaを読み出す。更に、入射強度算出部７は、反射率記憶部９から反射強度Iaに対応する既知の反射率Ｒaを読み出す。既知の反射率Ｒaは、図４に示すような一般的な分光反射特性のグラフから設定可能である。入射強度算出部７は、既知反射率領域における地盤反射強度Ia及び既知の反射率Ｒaを用いて、式（１）のように、図３に示すような、地盤に入射するレーザーの強度（以下、「入射強度」という。）I0を算出する。
【００２３】
I0＝Ia／Ra …（１）
例えば、地盤反射強度Iaの値を１２．２とし、使用したレーザーの波長が１．０６４［μm］程度であることを考慮して反射率Ｒaの値を０．１４とすると、入射強度I0の値は、８７．１２となる。算出された入射強度I0は、入射強度記憶部２７に格納される。
【００２４】
樹木頂点記憶部２３は、森林域Ａ内の樹木毎に、樹木頂点の座標を格納する（この樹木頂点の算出方法については後述する）。本発明の第１の実施の形態においては、一例として、図５に示すオルソフォト画像において、森林域Ａ内の樹木（樹冠）Ｔ１〜Ｔ６のＬＡＩを算出することとする。例えば、オルソフォト画像記憶部２０が、樹木（樹冠）Ｔ１〜Ｔ６の形状・座標等を、オルソフォト画像に対応させて格納していれば良い。なお、ＬＡＩの算出対象は樹木Ｔ１〜Ｔ６に限定されず、例えば森林域Ａ内のすべての樹木のＬＡＩを算出しても良い。樹木頂点記憶部２３は、図５に示した樹木Ｔ１〜Ｔ６毎に、樹木頂点Ｐi（ｉ＝１〜６）の座標（Ｘi，Ｙi）を格納する。
【００２５】
樹冠反射強度抽出部１３は、樹木頂点記憶部２３から樹木頂点Ｐiの座標（Ｘi，Ｙi）を読み出し、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）のラストパルスの強度Ifirstを、図３に示すような樹冠で反射したレーザーの反射強度（以下、「樹冠反射強度」という。）Ictとして抽出する。樹冠反射強度Ictは、樹冠反射強度記憶部２８に格納される。
【００２６】
樹冠通過強度抽出部１４は、樹木頂点記憶部２３に格納された樹木頂点Ｐiの座標（Ｘi，Ｙi）を参照して、反射強度データ記憶部２１から、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データの反射地点の高さ（Ｚ値）Ｚdを読み出す。更に、樹冠通過強度抽出部１４は、ＤＥＭ記憶部２４から、ＤＥＭの樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の地盤の高さＺhを読み出す。更に、樹冠通過強度抽出部１４は、反射地点の高さＺdと地盤の高さＺhとの差分値Ｚkを算出する。
【００２７】
更に、樹冠通過強度抽出部１４は、閾値記憶部２５から閾値（例えば、１ｍ）を読み出して、差分値Ｚkと比較する。差分値Ｚkが閾値以内の場合、レーザーが樹冠を通過し、地盤で反射したものと判断されるので、反射強度データのラストパルスの強度Ilastを、図３に示すような樹冠を通過して減衰した後に地盤で反射したレーザーの反射強度（以下、「樹冠通過強度」という。）Icbとして抽出する。図３に示すように、樹冠通過強度Icbは、樹冠を通過することで減衰しているので、樹冠内部の情報を有するといえる。抽出した樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。
【００２８】
地盤反射強度抽出部１１は、反射強度データ記憶部２１から、森林域Ａ内の反射強度データの反射地点の高さ（Ｚ値）Ｚdを読み出す。更に、地盤反射強度抽出部１１は、ＤＥＭ記憶部２４から、森林域Ａ内のＤＥＭの地盤の高さＺhを読み出す。更に、地盤反射強度抽出部１１は、反射地点の高さＺdと地盤の高さＺhとの差分値Ｚkを算出する。
【００２９】
更に、地盤反射強度抽出部１１は、閾値記憶部２５から閾値（例えば、１ｍ）を読み出して、差分値Ｚkと比較する。差分値Ｚkが閾値よりも小さい場合、レーザーが直接地盤に到達したものと判断されるので、反射強度データのファーストパルスの強度Ifirstを、図３に示すような地盤に直接達して反射したレーザーの反射強度（以下、「地盤反射強度」という。）Ibとして抽出する。抽出した樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。なお、樹冠通過強度Icbが森林域Ａ内の反射率Rbで得られることから、地盤反射強度Ibも森林域Ａから抽出している。
【００３０】
反射率算出部１２は、入射強度記憶部２７から入射強度I0を読み込み、地盤反射強度記憶部２６から森林域Ａ内の地盤反射強度Ibを読み込んで、式（２）を用いて、森林域Ａ内の反射率Ｒｂを算出する。
【００３１】
Rb＝Ib／I0 ・・・（２）
森林域Ａ内の反射率Ｒｂは、反射率記憶部２２に記憶される。
【００３２】
図６に示すように、樹木Ｔ１〜Ｔ６毎に、森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0が設定される。
【００３３】
葉面積指数算出部１５は、地盤反射強度記憶部２６、入射強度記憶部２７、及び樹冠反射強度記憶部２８、樹冠通過強度記憶部２９のそれぞれに格納された森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出して、樹木Ｔ１〜Ｔ６毎に、ＬＡＩを算出する。
【００３４】
地盤反射強度Ibは下の式で表される。
【００３５】
Ib＝I0・Rｂ …（３）
ベール−ランベルトの法則を参考にすると、葉間隙（gap fraction）は下の式（４）で表される。
【００３６】
Ｔ(ｚ)＝exp(−Ｇ・μ・ｚ) …（４）
ここで、Ｔは葉間隙、Ｇは葉の投影係数、μは葉群密度、ｚは梢端からの透過距離を表す。
【００３７】
樹冠通過強度Icbは、下の式（５）で表される。
【００３８】
Icb＝(I0−Ict) ・exp(-2G・μ・d)・Rb …（５）
ここで、ｄは樹冠厚さを表す。式（５）より、樹冠通過強度Icbと 森林域Ａ内の地盤反射強度Ibの比をとると式（６）のようになる。
【００３９】
Icb／Ib＝ ((I0−Ict)/ I0) ・ exp(-2G・μ・d) …（６）
葉面積指数は下の式（７）で表される。
【００４０】
ＬＡＩ＝μ・ｄ …（７）
葉を水平面に投影したと仮定して、Ｇ＝１とする。式（６）及び式（７）より、葉面積指数ＬＡＩは、式（８）で求めることができる。
【００４１】
ＬＡＩ＝-1/2 ・ ln( (Icb・I0) / (Ib・(I0 - Ict)) ) …（８）
葉面積指数算出部１５は、樹木Ｔ１〜Ｔ６毎に、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８及び樹冠通過強度記憶部２９から、森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出す。更に、葉面積指数算出部１５は、式（８）のように、ＬＡＩを算出する。算出したＬＡＩは、葉面積指数記憶部３１に格納される。葉面積指数記憶部３１は、図７に示すように、樹木Ｔ１〜Ｔ６毎に、ＬＡＩの値を格納する。出力部４は、ＬＡＩの値を、記録媒体に記録したり、画面に表示したり、印刷したりする。
【００４２】
次に、本発明の第１の実施の形態に係る葉面積指数算出方法の一例を、図８のフローチャートを参照しながら説明する。
【００４３】
（イ）ステップＳ１１において、地盤反射強度抽出部１１が、反射強度データ記憶部２１から、図３に示すような地盤に直接達して反射したレーザーの反射強度（地盤反射強度）Ibを抽出する。森林域Ａ内の地盤反射強度Ibは、地盤反射強度記憶部２６に格納される。
【００４４】
（ロ）ステップＳ１２において、地盤反射強度抽出部６が、反射強度データ記憶部２１から、図３に示すような既知反射率領域の地盤反射強度Iaを抽出し、地盤反射強度記憶部８に記憶させる。そして、入射強度算出部７が、地盤反射強度記憶部８から、既知反射率領域の地盤反射強度Iaを読み出し、更に反射率記憶部９から既知の反射率Ｒaを読み出して、式（１）を用いて、樹木Ｔ１〜Ｔ６の樹冠に到達する直前のレーザーの入射強度（入射強度）I0を算出する。入射強度I0は、入射強度記憶部２７に格納される。
【００４５】
（ハ）ステップＳ１３において、樹冠反射強度抽出部１３が、反射強度データ記憶部２１から、図３に示すような樹木Ｔ１〜Ｔ６の樹冠で反射したレーザーの反射強度（樹冠反射強度）Ictを、樹木Ｔ１〜Ｔ６毎に抽出する。樹冠反射強度Ictは、樹冠反射強度記憶部２８に格納される。
【００４６】
（ニ）ステップＳ１４において、樹冠通過強度抽出部１４は、反射強度データ記憶部２１から、図３に示すような、樹木Ｔ１〜Ｔ６の樹冠を通過した後に地盤で反射したレーザーの反射強度（樹冠通過強度）Icbを、樹木Ｔ１〜Ｔ６毎に抽出する。樹冠通過強度Icbは、樹冠通過強度記憶部２９に格納される。
【００４７】
（ホ）ステップＳ１５において、葉面積指数算出部１５が、地盤反射強度記憶部２６、入射強度記憶部２７、樹冠反射強度記憶部２８及び樹冠通過強度記憶部２９から森林域Ａ内の地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0を読み出して、式（７）を用いて、樹木Ｔ１〜Ｔ６毎にＬＡＩを算出する。ＬＡＩは、葉面積指数記憶部３１に格納される。
【００４８】
なお、ステップＳ１１〜Ｓ１４の手順は、それぞれ独立して行うことが可能であり、図８のフローチャートに示した順序に特に限定されない。
【００４９】
次に、図８のステップＳ１１の森林域Ａ内の地盤反射強度Ibの設定方法の詳細を、図９のフローチャートを参照しながら説明する。
【００５０】
（イ）ステップＳ２１において、表示部２が、オルソフォト画像を画面に表示している。利用者が、オルソフォト画像の森林域ＡをＬＡＩの算出対象とする指示を入力部３に入力する。表示処理部１０が、入力部３からの指示に応じて、森林域Ａを示す画面上のピクセル座標（Ｘｐ，Ｙｐ）を読み出す。そして、表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当するオルソフォト画像上の平面直角座標（Ｘo，Ｙo）を読み、地盤反射強度抽出部１１に引き渡す。
【００５１】
（ロ）ステップＳ２２において、地盤反射強度抽出部１１が、森林域Ａ内の反射強度データを引き当てる。ステップＳ２３において、地盤反射強度抽出部１１が、反射強度データに含まれるレーザーの反射地点の高さＺｄを抽出する。
【００５２】
（ハ）ステップＳ２４において、地盤反射強度抽出部１１が、ＤＥＭ記憶部２４から、森林域Ａ内のＤＥＭの地盤の高さＺｈを読み出す。ステップＳ２５において、地盤反射強度抽出部１１が、反射強度データの反射地点の高さＺｄと、ＤＥＭの地盤の高さＺｈとの差分値Ｚｋをとる。
【００５３】
（ニ）ステップＳ２６において、地盤反射強度抽出部１１が、差分値Ｚｋが閾値と比較する。差分値Ｚｋが閾値以上の場合は、ステップＳ２４に戻り、森林域Ａ内の他の反射強度データを引き当てる。他方、差分値Ｚｋが閾値以内場合は、ステップＳ２７に進む。
【００５４】
（ホ）ステップＳ２７において、地盤反射強度抽出部１１が、反射強度データのファーストパルスの強度Ifirstを、図３に示すような地盤反射強度Ibとして抽出する。ステップＳ２８において、地盤反射強度抽出部１１が、地盤反射強度Ibを地盤反射強度記憶部２６に記憶させる。
【００５５】
次に、図８のステップＳ１３の樹冠反射強度Ictの設定方法の詳細を、図１０のフローチャートを参照しながら説明する。
【００５６】
（イ）ステップＳ３１において、表示部２が、オルソフォト画像を表示している。表示処理部１０が、入力部３からの指示に応じて、表示部２の画面に表示された森林域Ａを示すピクセル座標（Ｘｐ，Ｙｐ）を読み出す。表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当するオルソフォト画像の平面直角座標（Ｘo，Ｙo）を読み、樹冠反射強度抽出部１３に引き渡す。
【００５７】
（ロ）ステップＳ３２において、樹冠反射強度抽出部１３が、平面直角座標（Ｘo，Ｙo）に基づいて、樹木頂点記憶部２３から森林域Ａ内に存在する樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）を読み出す。
【００５８】
（ハ）ステップＳ３３において、樹冠反射強度抽出部１３が、樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）のうち、樹木頂点Ｐ１に着目する。
【００５９】
（ニ）ステップＳ３４において、樹冠反射強度抽出部１３が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致、或いは最も近い座標（Ｘ，Ｙ）を有する反射強度データのファーストパルスの強度Ifirstを樹冠反射強度Ictとして抽出する。
【００６０】
（ホ）ステップＳ３５において、樹冠反射強度抽出部１３が、抽出した樹冠反射強度Ictを、樹木Ｔ１の樹木頂点Ｐ１に対応させて樹冠反射強度記憶部２８に記憶させる。
【００６１】
（ヘ）ステップＳ３６において、樹冠反射強度抽出部１３が、着目した樹木頂点Ｐｉが最後の樹木頂点Ｐ６であるか、即ち全ての樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictを設定したか判断する。全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictを抽出していなければ、ステップＳ３７においてｉを１インクリメントし、ステップＳ３４に戻り、次の樹木頂点Ｐ２に対して樹冠反射強度Ictを抽出する。この結果、全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）に対して樹冠反射強度Ictが抽出される。
【００６２】
次に、図８のステップＳ１４の樹冠通過強度Icbの設定方法の詳細を、図１１のフローチャートを参照しながら説明する。
【００６３】
（イ）ステップＳ４１において、表示部２が、オルソフォト画像を表示している。表示処理部１０が、入力部３からの指示に応じて、表示部２の画面に表示されたＬＡＩの算出対象となる樹木群の森林域Ａを示すピクセル座標（Ｘｐ，Ｙｐ）を読み出す。表示処理部１０が、ピクセル座標（Ｘｐ，Ｙｐ）に該当する平面直角座標（Ｘo，Ｙo）を読み、樹冠通過強度抽出部１４に引き渡す。
【００６４】
（ロ）ステップＳ４２において、樹冠通過強度抽出部１４が、樹木頂点記憶部２３から、森林域Ａ内に存在する樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）を読み出す。
【００６５】
（ハ）ステップＳ４３において、樹冠通過強度抽出部１４が、読み出した樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）のうち、樹木頂点Ｐ１に着目する。
【００６６】
（ニ）ステップＳ４４において、樹冠通過強度抽出部１４が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致、或いは最も近い座標（Ｘ，Ｙ）を有する反射強度データを引き当てる。
【００６７】
（ホ）ステップＳ４５において、樹冠通過強度抽出部１４が、反射強度データの反射地点の高さＺｄを抽出する。
【００６８】
（ヘ）ステップＳ４６において、樹冠通過強度抽出部１４が、樹木頂点Ｐ１の座標（Ｘ１，Ｙ１）に一致する、或いは最も近いＤＥＭの地盤の高さＺｈを読み出す。
【００６９】
（ト）ステップＳ４７において、樹冠通過強度抽出部１４が、反射強度データの反射地点の高さＺｄと、ＤＥＭの地盤の高さＺｈとの差分値Ｚｋをとる。
【００７０】
（チ）ステップＳ４８において、樹冠通過強度抽出部１４が、差分値Ｚｋが、閾値以上か判断する。閾値よりも大きい場合は、ステップＳ４９に進み樹木頂点Ｐｉの座標（Ｘｉ，Ｙｉ）に次に近い座標を有する反射強度データを引き当て、ステップS４４に戻る。他方、閾値以内の場合は、ステップＳ５０に進み、樹冠通過強度抽出部１４が、Ｚｄを有する反射強度データを引き当てる。
【００７１】
（リ）ステップＳ５１において、樹冠通過強度抽出部１４が、反射強度データのラストパルスの強度Ilastを樹冠通過強度Icbとして抽出する。
【００７２】
（ヌ）ステップＳ５２において、樹冠通過強度抽出部１４が、樹冠通過強度Icbを樹冠通過強度記憶部２９に樹木頂点Ｐ１と対応させて保存（設定）する。
【００７３】
（ル）ステップＳ５３において、樹冠通過強度抽出部１４が、着目したのが樹木頂点Ｐ６であるか、即ち全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbを設定したか判断する。全ての樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbを設定していなければ、ステップＳ５４に進み、次の樹木頂点Ｐ２を対象として、ステップＳ４４に戻る。この結果、全ての樹木Ｔ１〜Ｔ６の樹木頂点Ｐｉ（Ｐ１，Ｐ２，…Ｐ６）について樹冠通過強度Icbが抽出される。
【００７４】
本発明の第１の実施の形態によれば、従来、現地調査により局所的なＬＡＩを求めることしかできなかったが、航空レーザー計測でＬＡＩを推定できるので、現地に立ち入ることなく、広域かつ迅速にＬＡＩを推定することができる。
【００７５】
更に、地盤反射強度Ib、樹冠反射強度Ict、樹冠通過強度Icb及び入射強度I0等のレーザーの強度を活用し、樹冠内部の情報（レーザー反射強度の減衰量）を用いてＬＡＩを推定しているので、高精度なＬＡＩを推定することができる。
【００７６】
更に、樹木頂点記憶部２３に格納された樹木頂点Ｐiを用いて、樹冠反射強度抽出部１３が、樹木頂点Ｐiに一致、或いは最も近い座標（Ｘ，Ｙ）の樹冠反射強度Ictを抽出し、樹冠通過強度抽出部１４が、樹木頂点Ｐiに一致、或いは最も近い座標（Ｘ，Ｙ）の樹冠通過強度Icbを抽出することにより、ＬＡＩを推定することができる。樹木頂点Ｐiの決定方法は後述する。
【００７７】
このようにして求めたＬＡＩは森林のＣＯ2吸収量を推定する際のパラメータとして活用することもできる。
【００７８】
（第２の実施の形態）
本発明の第２の実施の形態として、樹木頂点の特定方法について説明する。図１３に示した樹木頂点記憶部２３に含まれる樹木頂点Ｐｉは、以下のように特定している。
【００７９】
本発明の第２の実施の形態は、樹冠を際立たせて再現し、単木の頂点を認識するものである。この単木の頂点の認識のために、新たに「尾根谷度」データを作成することによって樹冠形状を判別する。
【００８０】
尾根谷度は次の様に算出する。
【００８１】
本発明の第２の実施の形態では、ＤＥＭとＤＳＭの差のデータである樹木高のグリッドデータ（0.5m）を用いて、８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。検索範囲は1m・2m・5m・10mの4通り試算し、樹冠形状を適切に表現できた2mに決定した。そして、検索範囲2mにおける大まかな地形を表す、地上開度と地下開度とに挟まれる角度の2等分線とを尾根谷度（φ3）と定義する。
【００８２】
以下に図面を用いて樹木認識の実施の形態を説明する。
【００８３】
図１２は本発明の第２の実施の形態の樹木頂点認識装置の概略構成図である。図１２に示すように、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を記憶したデータベース１１０と、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model）を記憶したデータベース１１１と、特定地域のデジタル写真を記憶したデータベース１１２と、表示部１３と、樹木高算出部１１４と、樹木高データが記憶されるファイル（メモリ）１１５と、平滑フィルタ部１１６と（３ピクセル×３ピクセル）、尾根谷度算出部１１７と、尾根谷度が記憶されるファイル１１８と、表示処理部１１９と、局所最大フィルタ部１２０と、樹木頂点候補抽出部１２１と、樹木頂点候補データが記憶されるファイル１２２と、樹冠上部抽出部１２３と、樹冠上部データが記憶されるファイル１２４と、樹木頂点決定部１２５と、樹木頂点データが記憶されるファイル１２６と出力部１３０等を備えている。
【００８４】
データベース１１０は、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model：本発明の第２の実施の形態では０．５ｍ間隔の格子状のデータ）を記憶している。このグリッドデータは、地盤のＸ、Ｙ、Ｚ（地盤高Ｚａ）とを対応させて記憶させられている。つまり、レーザーデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて地面を復元する。
【００８５】
データベース１１１は、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model：本発明の第２の実施の形態では０．５ｍ間隔）を記憶している。このグリッドデータ（本発明の第２の実施の形態では０．５ｍ間隔の格子状のデータ）は、樹木表面のＸ、Ｙ、Ｚ（Ｚｂ：樹木表面高）とを対応させて記憶している。つまり、レーザーデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて表層を復元する。
【００８６】
すなわち、これらのグリッドデータは、格子間隔０．５ｍのＸ−Ｙ平面直角座標系に、緯度経度に対応するＸ値、Ｙ値とレーザー計測によって得られたＺ値を割り当てている。
【００８７】
樹木高算出部１１４は、データベース１１０のＤＳＭとデータベース１１１のＤＥＭとの差を求め、これをファイル１１５（メモリ）のＸ，Ｙ座標に割付ける。
【００８８】
平滑フィルタ部１１６は、Average filterとも呼ばれるものであり、ファイル１１５の樹木高グリッドデータＤｉに対して、各々のデータを中心とするＸ−Ｙ座標面３×３範囲（本発明の第２の実施の形態では１．５ｍ×１．５ｍ）で平均を求め、これをファイル１２６に保存する。
【００８９】
尾根谷度算出部１１７は、ファイル１２６に平均化された樹木高グリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２・・）が保存されると、樹木高グリッドデータＤｄｉが付加されている任意の平均化された樹木高グリッドデータＤｄｉ（以下樹木高グリッドデータＤｄｉという）から各々のデータを中心としたＸ−Ｙ座標面の任意半径の検索範囲（検索範囲は樹冠の立体形状の特徴を適切に表現できた2mに決定している。以下２ｍ範囲毎）で尾根谷度を求め、これをファイル１１８に保存する。
【００９０】
局所最大値フィルタ部１２０は、ファイル１１８の尾根谷度のデータＤｅｉに３×３範囲のフィルタを順次適用し、この結果をファイル１２７に保存する。
【００９１】
樹木頂点候補抽出部１２１は、ファイル１２７の樹冠立体形状の特徴データＤｆｉ（ＸＹ、尾根谷度）とファイル１１８のＤｅｉ（Ｘｐｉ、Ｙｐｉ、尾根谷度）とを比較し、同じものを頂点候補Ｄｇｉとしてファイル１２２に保存する。
【００９２】
樹冠上部抽出部１２３は、ファイル１１８の尾根谷度データ（樹冠立体形状の特徴データＤｅｉ）に森林の状態に応じた閾値を設けることによって樹冠の上部のデータＤｅｉを抽出し、これを樹冠上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル１２４に保存する。
【００９３】
樹木頂点決定部１２５は、ファイル１２４の樹冠の上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル１２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定し、これをファイル１２８に保存する。そして、樹冠上部に複数の仮の頂点が存在するときは、この樹冠上部で最も樹木高が高いものを単木の頂点として決定し、これをファイル１２９に保存する（Ｐｉ）。
【００９４】
表示処理部１１９は決定した頂点を例えば○印（△印、×印、番号等の識別）にして画面に表示する。このとき、写真画像と重ねて表示してもよい。
【００９５】
出力部１３０は、ファイル１２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。このとき、樹木高も表示してもよい。
【００９６】
（動作説明）
上記のように構成された樹木頂点特定装置について以下に動作を説明する。本発明の第２の実施の形態では表示部１３にデータベース１１０のＤＳＭ及びデータベース１１１のＤＥＭ並びにデジタル写真を重ね表示している。
【００９７】
図１３は本発明の第２の実施の形態の樹木頂点認識装置の概略動作を説明するフローチャートである。
【００９８】
本発明の第２の実施の形態では、樹木高算出部１１４が航空レーザー計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ７１）。また、地盤のグリッドデータ（以下、ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ７２）。
【００９９】
そして、樹木高算出部１１４がDSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）をとり樹木高グリッドデータＤｉとして、ファイル１１５に保存する（Ｓ７３）。
【０１００】
この樹木高のグリッドデータＤｉに対して、平滑フィルタ部１１６（Average filter16ともいう）が３×３範囲で順次平均化し、ファイル１２６に保存（樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・））する（Ｓ７４）。
【０１０１】
次に、尾根谷度算出部１１７が２ｍ範囲毎（Ｄｄｉ）に地上開度及び地下開度を算出する（Ｓ７５）。次に、地下開度及び地上開度を用いて尾根谷度算出部１１７が「尾根谷度」データを作成することによって樹冠の立体形状の特徴を判別できるようにした尾根谷画像を得てファイル１１８に保存する（Ｓ７６）。
【０１０２】
尾根谷度は次の様に算出する。樹木高のグリッドデータＤｄｉを用いて、８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。
【０１０３】
検索範囲は1m・2m・5m・10mの4通り試算し、樹冠形状を適切に表現できた2mに決定した。
【０１０４】
そして、検索範囲2mにおける大まかな地形を表す、地上開度と地下開度に挟まれる角度の2等分線を尾根谷度（φ3）と定義していく（図２０参照）。この尾根谷度のデータＤｅｉは表示処理部１１９によって表示（尾根谷度画像）される。
【０１０５】
さらに、尾根谷度のデータＤｅｉに局所最大値フィルタ部１２０を順次適用し（Ｓ７７）、これをファイル１２７に保存する。
【０１０６】
このとき、樹冠の大きさに限らず最小のエリアサイズである3×3範囲に固定して、樹頂点ではなく「樹頂点の候補点」を抽出し、ファイル１２２に樹木頂点候補Ｄｇｈ（ｄｇ１、・・）として保存する（Ｓ７８）。この３×３範囲に固定したのは、例えば５×５にすると、樹木頂点候補が存在しない樹冠上部が多くなるという実験結果を得たので、３×３範囲とした。
【０１０７】
一方で、「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった性質を利用して、樹木頂点候補抽出部１２１が尾根谷度データに森林の状態に応じた閾値を設けることによって樹冠の上部を抽出して（Ｓ７９）、冠上部抽出部１２４をラベリング化する（Ｓ８０）。
【０１０８】
さらに、樹木頂点決定部１２５が樹頂候補点と樹冠上部エリア毎に、樹木頂点候補点の中から最高樹木高を求め、その点を樹木頂点と決定（Ｓ８１）し、これを樹木頂点として抽出する（Ｓ８２）。この樹木頂点の決定については後述する。
【０１０９】
図１４は本発明の第２の実施の形態の樹木頂点認識装置の尾根谷度を求める処理を説明するフローチャートである。
【０１１０】
本発明の第２の実施の形態では、表示部１３にオルソフォト画像を表示している（Ｓ９０）。
【０１１１】
樹木高算出部１１４は、ＤＳＭ（Ｘ、Ｙ、Ｚ）及びＤＥＭ（Ｘ、Ｙ、Ｚ）とを読み込み（Ｓ９１）、差Ｚｃｉを求め（Ｓ９２）、これらを樹木高のグリッドデータＤｉ（Ｄ１、Ｄ２、・・・）としてファイル１１５に保存する。
【０１１２】
この樹木高データＤｉは、図１５に示すように、座標（Ｘ，Ｙ：平面直角座標）と樹木高Ｚｃｉとが対応させられたデータとなる。
【０１１３】
次に、平滑フィルタ部１１６がファイル１１５の多数の樹木高のグリッドデータＤｉに対してフィルタ１６で順次平均化（３×３範囲）し（Ｓ９３）、この樹木高のグリッドデータＤｄｉをファイル１２６に保存する。
【０１１４】
次に、尾根谷度算出部１１７は、個々のDdiにおいて２ｍ範囲で順次、尾根谷度を求めることによって樹冠の立体形状の特徴を判別する（Ｓ９４）。
【０１１５】
尾根谷度は、８方向に探索（木の種類によって１ｍ、５ｍ、１０ｍで行う場合もある）して地上開度と地下開度とを求め、地上開度と地下開度との差を２で割ることで求める（Ｓ９４）。つまり、樹木の冠部（葉、枝）の立体形状の特徴が求められる。
【０１１６】
ファイル１２６の樹木高のグリッドデータＤｄｉを用いて、図１６に示すように８方向（＝全方位）に探索して地上開度（φ1）と地下開度（φ2）を求める。検索範囲は樹冠の立体形状の特徴を適切に表現できた2mに決定している。
【０１１７】
図１７は任意の検索範囲における樹木高のグリッドデータＤｄｉを距離と高さからなる軸にプロットしたものである。また、図１７の黒点は０．５ｍ間隔のグリッドデータＤｄｉを示す（２ｍ範囲は図示せず）。
【０１１８】
つまり、図１７は着目する標本地点から任意半径の範囲で一方向を見た場合の大まかな樹冠形状を表している。そして、この図１７において、φ１が地上開度（地上角）であり、φ２が地下開度（地下角）である。そして、図１８（ａ）（ｂ）では地上開度と地下開度に挟まれる角度の2等分線と距離軸（水平軸）とで挟まれた角度を尾根谷度（φ3）と定義している。また、図１８の黒点は０．５ｍ間隔のグリッドデータＤｄｉを示す（２ｍ範囲は図示せず）
本形態では索範囲を２ｍとして、尾根谷度（φ3）を求め、尾根谷度のデータ（以下樹冠立体形状の特徴データＤｅｉ（Ｄｅｉ：Ｄｅ１、Ｄｅ２・・・）としてファイル１１８に保存する（Ｓ９６）。
【０１１９】
図１９が樹木高をY軸とし距離をＸ軸として軌跡Ａｉの樹冠形状を表現した説明図である。図２０が尾根谷度をY軸にし距離をＸ軸として軌跡Ａｉの樹冠の立体形状の特徴を表現した説明図である。
【０１２０】
つまり、図２０に示すように尾根谷度で軌跡Ａｉの樹冠の立体形状の特徴を表現すると、図１９より顕著に樹冠間の境界が現れる。図１９及び図２０の番号（１）〜（６）は写真から肉眼で読み取った樹冠エリアを示すものである。
【０１２１】
単木としての認識しやすさを向上するには、樹頂点とエッジ部分の値の差が大きいこと、樹頂点において値が同傾向を示すこと、エッジにおいて値が同傾向を示すこと等が挙げられる。図２０に示すように、従来までの「高さ」データより本手法の「尾根谷度」データの方が先にあげた条件を満たしていることがわかる。これにより、樹冠ポリゴンを自動作成する際の精度向上に寄与することができる。
【０１２２】
図２０においては、尾根谷度がマイナスの値となる点から急激に＋方向にあがり、そして急激にマイナスの値となる点までの範囲が単木の樹冠部であることが分かる。
【０１２３】
また、図２２は頂点候補と樹頂点の選定法を説明する図であり、各樹冠上部で樹木高が最も高い候補を選定することし示している。この樹木頂点の決定については詳細に後述する。
【０１２４】
次に、尾根谷度算出部１１７は、求めた尾根谷度を樹冠立体形状の特徴データＤｅｉとしてファイル１１８（ＸＹが対応させられている）に保存（図２１参照）する（Ｓ９６）。
【０１２５】
前述の地上開度は、各要素（Ｄｄｉ）に対し、その要素から8方向に探索し、考慮距離内で最大傾斜角を求める。8方向それぞれの最大傾斜角の平均を、天頂からの角度で表す。なお、地上開度の値は考慮範囲に依存し、判読を行いたい対象のスケールに合わせて設定する必要がある。
【０１２６】
また、地下開度は各要素に対し、その要素から8方向に探索し、考慮距離内で最小傾斜角を求める。8方向それぞれの最小傾斜角の平均を、天頂からの角度で表す。地上開度と同様に、値は考慮範囲に依存し、判読を行いたい対象のスケールに合わせて設定する必要がある。
【０１２７】
尾根谷度は、算出された地上開度から地下開度を引き、２で割ることによって計算する。尾根地形は値がプラスになり、谷地形は値がマイナスで表現される。
【０１２８】
そして、図１３の樹木頂点抽出処理に進む。樹木頂点抽出処理のステップＳ７でファイル１１８の樹冠立体形状の特徴データＤｅｉ（Ｄｅ１、Ｄｅ２・・）に対して３×３の局所最大値フィルタ部１２０（３×３範囲のウインドウを用いる）を適用して、これをファイル１２７に保存（Ｄｆｉ）する。
【０１２９】
ステップＳ７８で、樹木頂点候補抽出部１２１が樹冠エリア（(1)、(2)・・）に存在するピーク点の樹冠立体形状の特徴データＤｆｉ（図２２参照）を頂点候補（以下Ｄｇｉ）として抽出し、ファイル１２２に保存する（図２３参照）。具体的にはファイル１１８のデータの値とファイル１２７の同じ座標値のデータの値が同値の場合を「１」と、他を「０」としてファイル１２２に定義する（図２３参照）。
【０１３０】
ここで、局所最大フィルタ部１２０及び樹木頂点候補抽出部１２１の処理について説明を補充する。
【０１３１】
つまり、局所最大フィルタ部１２０（Local max filter２０ともいう）は、ファイル１１８の樹冠立体形状の特徴データＤｅｉ（Ｄｅ１、Ｄｅ２・・）に３×３範囲のウインドウをかけ、このウインドウ内で最大となる値をウィンドウの中央値に置き換える処理を順次行う。これらを樹冠立体形状の特徴データＤｆｉとしてファイル１２７に保存する。
【０１３２】
そして、樹木頂点候補算出部２１が樹冠立体形状の特徴データＤｆｉ（ＸＹ、尾根谷度）とファイル１１８のＤｅｉ（ＸＹ、尾根谷度）とを比較し、同じものを頂点候補Ｄｇｉとしてファイル１２２に保存している。
【０１３３】
一方、ステップＳ７９で樹冠上部抽出部１２３が「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった性質を利用して、ファイル１１８の尾根谷度データ（樹冠立体形状の特徴データＤｅｉ）に森林の状態に応じた閾値を設けることによって樹冠の上部のデータＤｅｉを抽出し、これを樹冠上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル１２４に保存する（図２４参照）。
【０１３４】
図２４（ａ）の樹冠の上部のデータＤｈｉは、閾値を超えるＤｅｉを抽出したものであり、例えば番号１のエリアのデータＤｅｉは、（Ｘｈ５、Ｙｈ５）を「１」、（Ｘｈ６、Ｙｈ６）を「１」・・（Ｘｈ１０、Ｙｈ１０）を「１」として示している。
【０１３５】
図２４の（ｂ）は図２４（ａ）を具体的に示すものであり、閾値を超えた頂点は「１」とする。「１」が連続している範囲が樹冠上部である。
【０１３６】
そして、ステップＳ８０でファイル１２４の中の頂点候補ＤｈｉのＸｈｉ、Ｙｈｉをポリゴン化（座標Ｘｈｉ、Ｙｈｉと色情報）して表示処理部１１９によってラベルリング（纏まったエリアを（上部）にラベルを振る）させる（図２５参照）。
【０１３７】
次に、ステップＳ８１で樹木頂点決定部１２５は、ファイル１２４の樹冠の上部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル１２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定する。そして、樹冠上部のエリア内に複数の仮の頂点が存在するときは、この樹冠上部のエリア（図１９、図２０においては、番号１、２・・のエリア）で最も樹木高が高いものを単木の頂点として決定し、これをファイル１２９に保存する（Ｐｉ）。表示処理部１１９は決定した頂点を画面に表示する（図２６参照）。すなわち、例えば、頂点候補Ｄｇｉは図２７に示すように表示されるが、頂点決定された画像は図２８に示すようになる。
【０１３８】
前述の頂点の決定について説明を補充する。例えば、図２９の（ａ）に示すように、ファイル１２４に第１列目が「００１１００１１１００」、第２列目が「０１１１１０１１１１０」、第３列目が「００１１０００１０００」、第４列目が「０００００１０００１１」という具合に頂点上部域が得られた場合に、図２９（ｂ）に示すように、「１」が隣接するまとまり（「０」が区切り）を１つの樹冠上部としてまとまり毎に通し番号を割り振っていく。
【０１３９】
そして、例えば、同一高の頂点が樹冠内に複数の頂点が存在することになると、一度、ファイル１２８のデータＤＱｉを表示させて、オペレータによって指定させ、指定されたものを頂点Ｐｉと決定してファイル１２９に保存する。
【０１４０】
次に、尾根谷度について説明を補充する。尾根谷度の特徴としては、図３０に示すようなsinカーブの地形を考えると、尾根地形は「尾根谷度＞0」、谷地形は「尾根谷度＜0」で表現されるとともに、中腹のところで「尾根谷度＝0」となることが挙げられる。
【０１４１】
ここで、わかりやすくするために、理想の樹木が図３０に示すような形状をしているとすると、樹頂点では「尾根谷度＞0」となっており、必ず「尾根谷度＝0」となる点を通過して樹冠のエッジに到達し、そこでは「尾根谷度＜0」となっていると考えられる。したがって、樹冠の上部は「尾根谷度＞0」で定義できると言える。
【０１４２】
このことを利用して、本業務では、図３１に示すように樹頂点を絞り込むこととした。つまり、尾根谷度の値を利用することによって樹冠上部を抽出し、この抽出した中に存在する樹頂候補点のうち、最も樹高（DSM-DEM）が高い候補点を樹頂点として選定することとした。
【０１４３】
尾根谷度の閾値について説明を補充する。
【０１４４】
尾根谷度に閾値を設定して樹冠上部を抽出する際、図３２に示すように閾値αの場合には樹木５が抽出できなくなってしまうことがわかる。一方で、閾値βの場合には樹木５と樹木６が一本として認識されてしまう。このように、閾値によって抽出される樹頂点が左右されしまう問題を抱えている。
【０１４５】
したがって、ここでは、樹冠上部を抽出するための最適な閾値について検討することとした。具体的には、尾根谷度の閾値を「10,15,20,25,30」の５種類に設定して、樹頂点の抽出結果が最も良いものを採用することとした。絞込み後の樹頂点の集計結果を表１に示す。
【表１】
【０１４６】
表１を見ると、「尾根谷度＞20」での結果が最もよかったことがわかる。これと比較して「尾根谷度＞10」での抽出本数が少ない理由は隣接する樹冠が一つとして認識されたためである。また、「尾根谷度＞30」での抽出本数が少ない理由は、閾値を高く設定したために抽出できなかった樹冠があったためである。
【０１４７】
従って、樹冠上部を抽出するための尾根谷度の閾値は「20」が最適であることがわかった。樹冠上部の抽出画像は図２５、絞込み後の樹頂点抽出画像は図２６に示している。
【０１４８】
上記のようにある山林の樹冠の頂点を算出することによって、単木を特定できるので山林の資産の管理等に適用できる。例えば、出力部１３０がファイル１２９の頂点のデータＰｉ（Ｘ、Ｙ、Ｚ）を出力して印刷させると、本数と高さが分かるので山林の価値等が自動的に把握できる。
【０１４９】
図３３においては、樹木頂点を算出した後で樹冠ポリゴンを作成している。樹冠ポリゴン作成は、ステップＳ８２で抽出された樹木頂点とファイル１１８の尾根谷度を用いてWetershedアリゴリズムにより、樹冠ポリゴン（図３４参照）を生成する（Ｓ８３、Ｓ８５）。
【０１５０】
また、ステップＳ８３の閾値の決定方法は下記の通りである。
【０１５１】
図１４の写真画像のようにオルソフォト画像において任意の線（軌跡Ａｉ）を引き、尾根谷度算出部が尾根谷度の横断図を作成し、これを表示する（図２０に該当）。
【０１５２】
次に、この横断線のスタートからの追加距離（図２０のX軸に相当）を追い、表示部のオルソ画像の樹冠と樹冠の境界を探す。
【０１５３】
そして、境界位置の追加距離を見て、その追加距離の尾根谷度（尾根谷度算出部が算出して表示させる）を図２０から作業者が読み取る。
【０１５４】
この行程を繰り返し、樹冠上部が形成される尾根谷度を読み取る。
【０１５５】
図２０の(1)では0と-50、(2)では-50と-20、(3)では-20と0、(4)では0と0、(5)では0と10、(6)では10と-10となった。
【０１５６】
この結果から、(1)〜(6)が分かれる尾根谷度10を閾値とする。
【０１５７】
閾値を0とすると(5)と(6)が一つの樹冠上部になるため、0は不採用とする。
【０１５８】
閾値を30とすると(2)では樹冠上部が2つ形成され、(5)が樹冠上部にならないので、これも不採用とする。よって、図２０では尾根谷度10程度が閾値としてふさわしいことになる。このようにして決定した閾値を用いている。
【０１５９】
なお、上記実施の形態では航空写真を用いて説明したが、樹木高のグリッドデータ又はＤＥＭ、ＤＳＭにグレイスケールを割り当ててもよい。
【０１６０】
例えば、地上開度と地下開度の差分画像をグレイに、傾斜を赤のチャンネルにいれて、擬似カラー画像を作成することにより、尾根や山頂部分を白っぽく、また谷や窪地を黒っぽく表現し、傾斜が急な部分ほど赤く表現する。このような表現の組み合わせにより、１枚でも立体感のある画像（立体赤色化マップともいう）を用いてもよい。
【０１６１】
例えば、地上開度の値が４０度から１２０度程度の範囲に収まる場合は、５０度から１１０度をグレイスケールに対応させ、２５５諧調に割り当てる。つまり、尾根の部分（凸部）の部分ほど地上開度の値が大きいので、色が白くなる。
【０１６２】
また、本発明の第２の実施の形態は、図３５に示すようなコンピュータシステムVPSによって実現するのが好ましい。ＶＰＳ１は、ワークステーション、プロセッサ、マイクロコンピュータ、ロジック、レジスタ等の適宜な組み合わせからなる中央情報処理装置（CPU）と、この中央情報処理装置に必要な制御・操作情報を入力するキーボード（KB）、マウス、対話型ソフトスイッチ、外部通信チャネル等を含む情報入力部と、中央情報処理装置から出力された情報を広義な意味で表示・伝送するディスプレイ、プリンタ、外部通信チャネル等を含む情報出力部と、中央情報処理装置に読み込まれるオペレーティングシステム、アプリケーションプログラム等の情報が格納されたロム（ＲＯＭ）と、中央情報処理装置で随時処理すべき情報及び中央情報処理装置から随時書き込まれる情報を格納するラム（ＲＡＭ）等とを備える。ＲＯＭ、ＲＡＭを適宜統合、細分化することは差し支えない。
【０１６３】
前述のアプリケーションプログラムは、本発明の第２の実施の形態の、樹木高算出部１１４と、尾根谷度算出部１１７と、表示処理部１１９と、局所最大フィルタ部１２０（Ｌｏｃａｌ max ｆｉｌｔｅｒともいう）と、樹木頂点候補抽出部１２１と、樹冠上部抽出部１２３と、樹木頂点決定部１２５と、出力部１３０等である。
【０１６４】
本実施形態では、開度という概念を用いている。開度は当該地点が周囲に比べて地上に突き出ている程度及び地下に食い込んでいる程度を数量化したものである。つまり、地上開度は、着目する標本地点から距離Ｌの範囲内で見える空の広さを表しており、また地下開度は逆立ちをして地中を見渡す時、距離Ｌの範囲における地下の広さを表している。
【０１６５】
開度は距離Ｌと周辺地形に依存している。一般に地上開度は周囲から高く突き出ている地点ほど大きくなり、山頂や尾根では大きな値をとり窪地や谷底では小さい。逆に地下開度は地下に低く食い込んでいる地点ほど大きくなり、窪地や谷底では大きな値をとり山頂や尾根では小さい。
【０１６６】
また、開度図は計算距離の指定によって、地形規模に適合した情報抽出が可能であり、方向性及び局所ノイズに依存しない表示が可能である。
【０１６７】
つまり、尾根線及び谷線の抽出に優れており、豊富な地形・地質情報が判読できるものであり、一定範囲のＤＥＭデータ上において、設定した当該地点Ａから８方向のいずれか一方を見たときに最大頂点となる点Ｂを結ぶ直線Ｌ１と、水平線とがなす角度ベクトルθｉを求める。この角度ベクトルの求め方を８方向に渡って実施し、これらを平均化したものを地上開度と称し、一定範囲のＤＥＭデータ上（地表面：立体）に空気層を押し当てた立体を裏返した反転ＤＥＭデータの当該地点Ａから８方向のいずれか一方を見たときに最大頂点となる点Ｃ（一番深い所に相当する）を結ぶ直線Ｌ２と、水平線とがなす角度を求める。この角度を８方向に渡って求めて平均化したのを地下開度ψｉと称している。
【０１６８】
すなわち、地上開度は、着目点から一定距離までの範囲に含まれるＤＥＭデータ上において、８方向毎に地形断面を生成し、それぞれの地点と着目点を結ぶ線Ｌ１の傾斜の最大値（鉛直方向から見たとき）を求める。このような処理を８方向に対して行う。傾斜の角度は天頂からの角度（平坦なら９０度、尾根や山頂では９０度以上、谷底や窪地では９０度以下） また、地下開度は、反転ＤＥＭデータの着目点から一定距離までの範囲において、８方向毎に地形断面を生成し、それぞれの地点と着目点を結ぶ線の傾斜の最大値（地表面の立体図において鉛直方向からＬ２を見たときには最小値）を求める。このような処理を８方向に対して行う。
【０１６９】
地表面の立体図において鉛直方向からＬ２を見たときの角度ψｉは、平坦なら９０度、尾根や山頂では９０度以下、谷底や窪地では９０度以上である。
【０１７０】
つまり、地上開度と地下開度は、２つの基本地点Ａ（ｉA，ｊA，ＨA）とＢ（ｉB，ｊB，ＨB）を考える。
【０１７１】
ＡとＢの距離は
P = {(iA − iB)2 + (jA − jB)2}1/2 …(1)
となる。
【０１７２】
標高０ｍを基準として、標本地点のＡとＢの関係を示したものである。標本地点Ａから標本地点Ｂを見た仰角θは
θ＝ｔａｎ-1｛（ＨB−ＨA）／Ｐ｝
で与えられる。θの符号は(1)ＨA＜ＨB の場合には正となり、(2)ＨA＞ＨB の場合には負となる。
【０１７３】
着目する標本地点から方位Ｄ距離Ｌの範囲内にある標本地点の集合を DＳL と記述して、これを「着目する標本地点のＤ−Ｌ集合」を呼ぶことにする。ここで、
DβL ：着目する標本地点の DＳL の各要素に対する仰角のうちの最大値
DδL ：着目する標本地点の DＳL の各要素に対する仰角のうちの最小値
として、次の定義をおこなう。
【０１７４】
定義１：着目する標本地点のＤ−Ｌ集合の地上角及び地下角とは、各々
DφL ＝９０− DβL
及び
DψL ＝９０＋ DδL
を意味するものとする。
【０１７５】
DφL は着目する標本地点から距離Ｌ以内で方位Ｄの空を見ることができる天頂角の最大値を意味している。一般に言われる地平線角とはＬを無限大にした場合の地上角に相当している。また、DψL は着目する標本地点から距離Ｌ以内で方位Ｄの地中を見ることができる天底角の最大値を意味している。Ｌを増大させると、 DＳL に属する標本地点の数は増加することから、 DβL に対して非減少特性を持ち、逆に DδL は非増加特性を持つ。したがって DφL 及びDψLは共にＬに対して非増加特性を持つことになる。
【０１７６】
測量学における高角度とは、着目する標本地点を通過する水平面を基準にして定義される概念であり、θとは厳密には一致しない。また地上角及び地下角を厳密に議論しようとすれば、地球の曲率も考慮しなければならず、定義１は必ずしも正確な記述ではない。定義１はあくまでもＤＥＭを用いて地形解析をおこなうことを前提として定義された概念である。
【０１７７】
地上角及び地下角は指定された方位Ｄについての概念であったが、これを拡張したものとして、次の定義を導入する。
【０１７８】
定義II：着目する標本地点の距離Ｌの地上開度及び地下開度とは、各々 ΦL＝（0φL ＋45φL ＋90φL ＋135φL ＋180φL ＋225φL ＋270φL ＋315φL ）／８
及び
ΨL＝（０ψL ＋45ψL ＋90ψL ＋135ψL ＋180ψL ＋225ψL ＋270ψL ＋315ψL ）／８
を意味するものとする。
【０１７９】
地上開度は着目する標本地点から距離Ｌの範囲内で見える空の広さを表しており、また地下開度は逆立ちをして地中を見渡す時、距離Ｌの範囲における地下の広さを表している。
【０１８０】
（その他の実施の形態）
上記のように、本発明は第１及び第２の実施の形態によって記載したが、この開示の一部をなす論述及び図面はこの発明を限定するものであると理解すべきではない。この開示から当業者には様々な代替実施の形態、実施例及び運用技術が明らかとなろう。
【０１８１】
例えば、本発明の第１の実施の形態では、樹冠反射強度抽出部１３が、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ifirstを、樹冠反射強度Ictとして抽出する例を説明したが、樹木頂点Ｐiの代わりに、本発明の第２の実施の形態において説明した樹冠上部データ抽出部１２３により抽出された樹冠上部のデータＤhiを用いても良い。即ち、樹冠反射強度抽出部１３が、樹冠上部のデータＤhiを用いて、樹冠上部付近の座標（Ｘ，Ｙ）のラストパルスの強度Ifirstを、樹冠反射強度Ictとして抽出しても良い。
【０１８２】
また、本発明の第１の実施の形態では、樹冠通過強度抽出部１４が、樹木頂点Ｐiの座標（Ｘi，Ｙi）に一致、或いは最も近い座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ilastを、樹冠通過強度Icbとして抽出する例を説明したが、樹木頂点Ｐiの代わりに、本発明の第２の実施の形態において説明した樹冠上部データ抽出部１２３により抽出された樹冠上部のデータＤhiを用いても良い。即ち、冠通過強度抽出部１４が、樹冠上部のデータＤhiを用いて、樹冠上部付近の座標（Ｘ，Ｙ）の反射強度データのラストパルスの強度Ilastを、樹冠通過強度Icbとして抽出しても良い。
【０１８３】
また、本発明の第２の実施の形態において樹木頂点の決定方法を説明したが、本発明の第１の実施の形態において用いる樹木頂点は、本発明の第２の実施の形態において説明した決定方法によるものでなくても良く、決定方法は特に限定されない。例えば、樹木頂点の特定方法としては、写真画像等の樹木の頂点を目視により特定する手法や、航空レーザー計測により取得されたデータを用いてＤＳＭとＤＥＭを作成し、ＤＳＭとＤＥＭの差分をとることによって樹木高のグリッドデータを算出し、樹木高や樹木高の極大値を判別指標として樹木頂点を特定する手法を用いても良い。
【０１８４】
このように、本発明はここでは記載していない様々な実施の形態等を含むことは勿論である。したがって、本発明の技術的範囲は上記の説明から妥当な特許請求の範囲に係る発明特定事項によってのみ定められるものである。
【図面の簡単な説明】
【０１８５】
【図１】本発明の第１の実施の形態に係る具体的なハードウエア構成図である。
【図２】本発明の第１の実施の形態に係る葉面積指数算出装置の概略構成図である。
【図３】レーザー反射強度を説明する説明図である。
【図４】分光特性（反射率）のグラフである。
【図５】算出エリアを説明するための平面図である。
【図６】抽出した樹冠内部情報（レーザー反射強度）を示す表である。
【図７】算出した葉面積指数を示す表である。
【図８】本発明の第１の実施の形態に係る葉面積指数算出方法のフローチャートである。
【図９】地盤反射強度の設定方法のフローチャートである。
【図１０】樹冠反射強度の設定方法のフローチャートである。
【図１１】樹冠通過強度の設定方法のフローチャートである。
【図１２】本発明の第２の実施の形態の樹木頂点認識装置の概略構成図である。
【図１３】本発明の第２の実施の形態の樹木頂点認識装置の概略動作を説明するフローチャートである。
【図１４】本発明の第２の実施の形態の樹木頂点認識装置の尾根谷度を求める処理を説明するフローチャートである。
【図１５】樹木高データＤiの説明図である。
【図１６】尾根谷度の算出を説明する説明図である。
【図１７】尾根谷度の算出を説明する説明図である。
【図１８】尾根谷度の算出を説明する説明図である。
【図１９】樹木高をＹ軸とし距離をＸ軸として軌跡Ａｉの樹冠形状を表現した説明図である。
【図２０】尾根谷度をＹ軸にし距離をＸ軸として軌跡Ａｉの樹冠の立体形状の特徴を表現した説明図である。
【図２１】樹木高グリッドデータの説明図である。
【図２２】尾根谷度による頂点候補の算出を説明する説明図である。
【図２３】頂点候補Ｄgiの説明図である。
【図２４】樹冠の上部のデータＤhiの説明図である。
【図２５】樹冠上部のラベリングを説明する説明図である。
【図２６】決定した頂点を画面の説明図である。
【図２７】頂点決定までの説明図である。
【図２８】頂点決定までの説明図である。
【図２９】頂点決定までの説明図である。
【図３０】尾根谷度の説明図である。
【図３１】尾根谷度の説明図である。
【図３２】閾値と樹冠部の関係を説明する説明図である。
【図３３】他の実施の形態を説明するフローチャートである。
【図３４】他の実施の形態によって得られた樹冠のポリゴンを説明する説明図である。
【図３５】本発明の第２の実施の形態の具体的なハードウエア構成図である。
【符号の説明】
【０１８６】
１…ＣＰＵ
２…表示部
３…入力部
４…出力部
５…記憶装置
１０…表示処理部
１１…地盤反射強度抽出部
１２…I0算出部
１３…樹冠反射強度抽出部
１４…樹冠通過強度抽出部
１５…葉面積指数算出部
１６…平滑フィルタ（アベレージフィルタ）
２０…オルソフォト画像記憶部
２１…反射強度データ記憶部
２２…反射率記憶部
２３…樹木頂点記憶部
２４…ＤＥＭ記憶部
２５…閾値記憶部
２６…地盤反射強度記憶部
２７…入射強度記憶部
２８…樹冠反射強度記憶部
２９…樹冠通過強度記憶部
３１…葉面積指数記憶部
１１０…データベース
１１１…データベース
１１２…データベース
１１３…表示部
１１４…樹木高算出部
１１６…平滑フィルタ
１１７…尾根谷度算出部
１１８…ファイル
１１９…表示処理部
１２０…局所最大値フィルタ
１２１…樹木頂点候補抽出部
１２２…ファイル
１２３…樹冠上部データ抽出部
１２４…ファイル
１２５…樹木頂点決定部
１２６…ファイル
１２７…ファイル
１２８…ファイル
１２９…ファイル
１３０…出力部
【特許請求の範囲】
【請求項１】
上空から地上へレーザーを照射し、前記地上で反射した前記レーザーの反射強度データを記憶するステップと、
地盤の反射率を記憶するステップと、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出するステップと、
記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出するステップと、
記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出するステップと、
記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出するステップと、
前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出するステップ
とを含むことを特徴とする葉面積指数算出方法。
【請求項２】
前記樹冠に到達する直前のレーザーの入射強度を算出するステップが、
記憶した前記反射強度データから、前記地盤で反射したレーザーの反射強度を抽出し、
前記地盤で反射したレーザーの反射強度を前記地盤の反射率で割って、前記樹冠に到達する直前のレーザーの入射強度を算出する
ことを特徴とする請求項１に記載の葉面積指数算出方法。
【請求項３】
前記反射強度データが、前記レーザーのファーストパルスの強度、前記レーザーのラストパルスの強度及び前記レーザーの反射地点の高さを含み、
前記地盤の高さを記憶するステップと、
特定の閾値を記憶するステップと、
前記反射地点の高さと前記地盤の高さとの差分値と前記閾値を比較するステップと、
前記比較の結果、前記差分値が前記閾値以内である前記反射強度データに含まれる前記ラストパルスの強度を、前記樹冠を通過したレーザーの反射強度として抽出するステップと、
前記比較の結果、前記差分値が前記閾値よりも小さい前記反射強度データに含まれる前記ファーストパルスの強度を、前記地盤で反射したレーザーの反射強度として抽出するステップ
とを含むことを特徴とする請求項１又は２に記載の葉面積指数算出方法。
【請求項４】
前記地盤で反射したレーザーの反射強度を、前記樹木の樹冠に到達する直前のレーザーの入射強度で割ることにより、前記反射率を算出することを特徴とする請求項１〜３のいずれか１項に記載の葉面積指数算出方法。
【請求項５】
（プログラム。請求項５〜８は請求項１〜４と同じ内容です。）
コンピュータを、
上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手段と、
地盤の反射率を記憶する手段と、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手段と、
記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手段と、
記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出する手段と、
記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出する手段と、
前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出する手段
として機能させるためのプログラム。
【請求項６】
前記コンピュータを、
記憶した前記反射強度データから、前記地盤で反射したレーザーの反射強度を抽出する手段と、
前記地盤で反射したレーザーの反射強度を前記地盤の反射率で割って、前記樹冠に到達する直前のレーザーの入射強度を算出する手段
として機能させるための請求項５に記載のプログラム。
【請求項７】
前記反射強度データが、前記レーザーのファーストパルスの強度、前記レーザーのラストパルスの強度及び前記レーザーの反射地点の高さを含み、
前記コンピュータを、
前記地盤の高さを記憶する手段と、
特定の閾値を記憶する手段と、
前記反射地点の高さと前記地盤の高さとの差分値と前記閾値を比較する手段と、
前記比較の結果、前記差分値が前記閾値以上である前記反射強度データに含まれる前記ラストパルスの強度を、前記樹冠を通過したレーザーの反射強度として抽出する手段と、
前記比較の結果、前記差分値が前記閾値よりも小さい前記反射強度データに含まれる前記ファーストパルスの強度を、前記地盤で反射したレーザーの反射強度として抽出する手段
として機能させるための請求項５又は６に記載のプログラム。
【請求項８】
前記コンピュータを、
前記地盤で反射したレーザーの反射強度を、前記樹木の樹冠に到達する直前のレーザーの入射強度で割ることにより、前記反射率を算出する手段
として機能させるための請求項５〜７のいずれか１項に記載のプログラム。
【請求項９】
上空から地上へレーザーを照射するレーザー照射手段と、
前記地上で反射した前記レーザーの反射強度データを記憶する反射強度記憶手段と、
地盤の反射率を記憶する反射率記憶手段と、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出し、記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出し、記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出し、記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出し、前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出する演算手段と、
前記葉面積指数を出力する出力手段
とを備える葉面積指数算出装置。
【請求項１】
上空から地上へレーザーを照射し、前記地上で反射した前記レーザーの反射強度データを記憶するステップと、
地盤の反射率を記憶するステップと、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出するステップと、
記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出するステップと、
記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出するステップと、
記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出するステップと、
前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出するステップ
とを含むことを特徴とする葉面積指数算出方法。
【請求項２】
前記樹冠に到達する直前のレーザーの入射強度を算出するステップが、
記憶した前記反射強度データから、前記地盤で反射したレーザーの反射強度を抽出し、
前記地盤で反射したレーザーの反射強度を前記地盤の反射率で割って、前記樹冠に到達する直前のレーザーの入射強度を算出する
ことを特徴とする請求項１に記載の葉面積指数算出方法。
【請求項３】
前記反射強度データが、前記レーザーのファーストパルスの強度、前記レーザーのラストパルスの強度及び前記レーザーの反射地点の高さを含み、
前記地盤の高さを記憶するステップと、
特定の閾値を記憶するステップと、
前記反射地点の高さと前記地盤の高さとの差分値と前記閾値を比較するステップと、
前記比較の結果、前記差分値が前記閾値以内である前記反射強度データに含まれる前記ラストパルスの強度を、前記樹冠を通過したレーザーの反射強度として抽出するステップと、
前記比較の結果、前記差分値が前記閾値よりも小さい前記反射強度データに含まれる前記ファーストパルスの強度を、前記地盤で反射したレーザーの反射強度として抽出するステップ
とを含むことを特徴とする請求項１又は２に記載の葉面積指数算出方法。
【請求項４】
前記地盤で反射したレーザーの反射強度を、前記樹木の樹冠に到達する直前のレーザーの入射強度で割ることにより、前記反射率を算出することを特徴とする請求項１〜３のいずれか１項に記載の葉面積指数算出方法。
【請求項５】
（プログラム。請求項５〜８は請求項１〜４と同じ内容です。）
コンピュータを、
上空から地上へレーザーを照射し、地上で反射したレーザーの反射強度データを記憶する手段と、
地盤の反射率を記憶する手段と、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出する手段と、
記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出する手段と、
記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出する手段と、
記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出する手段と、
前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出する手段
として機能させるためのプログラム。
【請求項６】
前記コンピュータを、
記憶した前記反射強度データから、前記地盤で反射したレーザーの反射強度を抽出する手段と、
前記地盤で反射したレーザーの反射強度を前記地盤の反射率で割って、前記樹冠に到達する直前のレーザーの入射強度を算出する手段
として機能させるための請求項５に記載のプログラム。
【請求項７】
前記反射強度データが、前記レーザーのファーストパルスの強度、前記レーザーのラストパルスの強度及び前記レーザーの反射地点の高さを含み、
前記コンピュータを、
前記地盤の高さを記憶する手段と、
特定の閾値を記憶する手段と、
前記反射地点の高さと前記地盤の高さとの差分値と前記閾値を比較する手段と、
前記比較の結果、前記差分値が前記閾値以上である前記反射強度データに含まれる前記ラストパルスの強度を、前記樹冠を通過したレーザーの反射強度として抽出する手段と、
前記比較の結果、前記差分値が前記閾値よりも小さい前記反射強度データに含まれる前記ファーストパルスの強度を、前記地盤で反射したレーザーの反射強度として抽出する手段
として機能させるための請求項５又は６に記載のプログラム。
【請求項８】
前記コンピュータを、
前記地盤で反射したレーザーの反射強度を、前記樹木の樹冠に到達する直前のレーザーの入射強度で割ることにより、前記反射率を算出する手段
として機能させるための請求項５〜７のいずれか１項に記載のプログラム。
【請求項９】
上空から地上へレーザーを照射するレーザー照射手段と、
前記地上で反射した前記レーザーの反射強度データを記憶する反射強度記憶手段と、
地盤の反射率を記憶する反射率記憶手段と、
記憶した前記反射強度データから、樹木の樹冠で反射したレーザーの反射強度を抽出し、記憶した前記反射強度データから、前記樹冠を通過して地盤で反射したレーザーの反射強度を抽出し、記憶した前記反射強度データ及び反射率に基づいて、前記樹冠に到達する直前のレーザーの入射強度を算出し、記憶した前記反射強度データから、地盤で反射したレーザーの反射強度を抽出し、前記樹冠に到達する直前のレーザーの入射強度、前記樹冠を通過したレーザーの反射強度、前記樹冠で反射したレーザーの反射強度、及び前記地盤で反射したレーザーの反射強度を用いて、前記樹木の葉面積指数を算出する演算手段と、
前記葉面積指数を出力する出力手段
とを備える葉面積指数算出装置。
【図１】
【図２】
【図３】
【図４】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３５】
【図５】
【図１４】
【図２５】
【図２６】
【図２７】
【図２８】
【図３４】
【図２】
【図３】
【図４】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３５】
【図５】
【図１４】
【図２５】
【図２６】
【図２７】
【図２８】
【図３４】
【公開番号】特開２００８−１１１７２５（Ｐ２００８−１１１７２５Ａ）
【公開日】平成２０年５月１５日（２００８．５．１５）
【国際特許分類】
【出願番号】特願２００６−２９４８２０（Ｐ２００６−２９４８２０）
【出願日】平成１８年１０月３０日（２００６．１０．３０）
【出願人】（０００００３６８７）東京電力株式会社 (2,580)
【出願人】（５９１０７４１６１）アジア航測株式会社 (48)
【Ｆターム（参考）】
【公開日】平成２０年５月１５日（２００８．５．１５）
【国際特許分類】
【出願日】平成１８年１０月３０日（２００６．１０．３０）
【出願人】（０００００３６８７）東京電力株式会社 (2,580)
【出願人】（５９１０７４１６１）アジア航測株式会社 (48)
【Ｆターム（参考）】
[ Back to top ]