樹木頂点認識方法及び樹木頂点認識装置並びに樹木頂点認識のプログラム
【課題】自動的に単木の頂点を認識する方法および認識装置を提供する。
【解決手段】樹木高算出部１４、検索範囲決定部１７、第２の樹冠形状指数算出部４１、第２の樹木頂点決定部４３、閾値算出部３５、局所領域設定部３６等を備え、検索範囲に対応する樹木高グリッドデータをファイル１５から読み込み、このグリッドデータにおけるアスペクト比の最大、最小値を基に樹冠形状をより明瞭にする所定の値に変更し、この変更された最大値、最小値を用い新たな地上開度、地下開度を求め、これらの開度から算出される樹冠形状指数を平均化し、平均化された樹冠形状指数を局所領域又はガウシアンフィルタで平滑化し、検索範囲における樹冠形状指数として算出する。そして、予め設定された局所領域の画像の平均値と標準偏差との和を閾値Ｈｇｉとして求め、閾値Ｈｇｉより大きい値の樹冠形状指数を樹冠部として抽出し、この樹冠部から最適な単木の頂点を決定する。
【解決手段】樹木高算出部１４、検索範囲決定部１７、第２の樹冠形状指数算出部４１、第２の樹木頂点決定部４３、閾値算出部３５、局所領域設定部３６等を備え、検索範囲に対応する樹木高グリッドデータをファイル１５から読み込み、このグリッドデータにおけるアスペクト比の最大、最小値を基に樹冠形状をより明瞭にする所定の値に変更し、この変更された最大値、最小値を用い新たな地上開度、地下開度を求め、これらの開度から算出される樹冠形状指数を平均化し、平均化された樹冠形状指数を局所領域又はガウシアンフィルタで平滑化し、検索範囲における樹冠形状指数として算出する。そして、予め設定された局所領域の画像の平均値と標準偏差との和を閾値Ｈｇｉとして求め、閾値Ｈｇｉより大きい値の樹冠形状指数を樹冠部として抽出し、この樹冠部から最適な単木の頂点を決定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、レーザ計測で得た多数のデータ（ＤＳＭ、ＤＥＭ）から単木の頂点を認識する樹木頂点認識装置に関する。
【背景技術】
【０００２】
森林のCO2吸収、水源涵養、木材生産などの様々な機能を向上させるためには森林管理が必要である。森林管理を行うに当たっては、樹木本数などの林分パラメータが必要となる。
【０００３】
例えば、山地の特定エリアの唐松の本数を特定する場合には、写真画像等のそのエリア内の唐松の樹木の頂点を目視で特定し、これらの総数を求めていた。
【０００４】
また、レーザ測距儀（以下、航空レーザ）を活用して樹木本数を推定する場合もある。
【０００５】
具体的には、航空レーザ計測により取得されたデータを用いて表層のグリッドデータ（以下、ＤＳＭ：Digital Surface Model）と地盤のグリッドデータ（以下、ＤＥＭ：Digital Elevation Model）を作成する。さらに、DSMとDEMデータの差分をとることによって樹木高のグリッドデータを算出し、樹木高や樹木高の極大値を判別指標として単木抽出を行っていた。
【０００６】
また、特開２００４−３３１４９号公報（特許文献１）は、施業林分における立木の直径、本数及び配置の調査を自動的に行うものである。
【特許文献１】特開２００４−３３１４９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、写真画像から目視によって樹木本数を求めた場合は、非常に時間がかかる。特に、エリアが広くなるとそれだけ時間がかかる。また、写真画像から樹木を１本毎に目視で正しく特定していくには、経験が必要である。
【０００８】
また、従来のＤＥＭ、ＤＳＭによって単木抽出を行う手法は、樹冠形状が複雑であることから、一つの樹冠内には複数の極大値が存在しやすく、従来の手法で単木抽出処理を行うと実本数に比べ過剰に抽出されるデメリットがある。
【０００９】
本発明は以上の課題を解決するためになされたもので、自動的に単木の頂点を正しく認識する樹木認識装置を得ることを目的とする。
【課題を解決するための手段】
【００１０】
本発明の樹木頂点認識方法は、上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識方法である。
【００１１】
コンピュータが、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶手段に記憶するステップと、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出ステップと、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出ステップと、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出ステップと、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識ステップと
を行うことを要旨とする。
【発明の効果】
【００１２】
以上のように本発明によれば、樹冠形状指数によって、自動的に精度よく単木の頂点を認識できる。
【００１３】
このため、広い森林等を人手によって樹木一本毎に検索しなくともよいので、例えば森林の財産管理等に適用できる。
【発明を実施するための最良の形態】
【００１４】
＜実施の形態１＞
初めに経緯を説明する。
【００１５】
尾根谷度の算出は、次の様に行う。樹木高のグリッドデータを用いて、設定した対象範囲内（検索範囲）における天頂から地表（樹冠含む）に接するまでの角度と鉛直下方から地表に接するまでの角度を求める。そして、これらの角度を例えば、8方向で計算し、それぞれ平均する。その結果をそれぞれ地上開度（φ1）と地下開度（φ2）とする。地上開度と地下開度にはさまれる角度の2 等分線と水平線との角度を尾根谷度（φ3）とする（図１参照）。検索範囲は1m・2m・5m・10mの4 通り試算し、樹冠形状を適切に表現できた2m に決定した。
【００１６】
尾根谷度は大まかな樹冠形状を表現するものであり、樹高に比べ、樹木頂点と樹冠端部の値の差が大きく、樹木頂点において値が同じような傾向を示す（図２参照）。これにより画一的な閾値を用いて、樹冠部域の抽出が可能で、ここに尾根谷度の有用性があると考える。
【００１７】
次に、尾根谷度に局所最大値フィルタを適用し、局所領域内の最大値を中央画素値に置き換える。この際、局所領域を3×3 に設定した。このフィルタ処理後の尾根谷度と処理前の尾根谷度を比較し、処理前後で尾根谷度が同値であったときに、その地点を「樹頂候補点」として抽出する。
【００１８】
一方で、「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった尾根谷度の性質を利用して、尾根谷度データに森林の状態に応じた閾値を設けることによって樹冠部を抽出する。そして、樹頂候補点と樹冠部領域を重ね合わせて、各樹冠部領域内において樹高が最も高い樹頂候補点を樹木頂点として抽出する。
【００１９】
しかしながら、抽出された樹木頂点の総数を現地の立木本数として検証していたが、抽出された樹木頂点は実際の立木本数よりも過小になる場合があることがわかった。すなわち、これまでの手法では抽出本数が過小になる可能性がある。
【００２０】
例えば、図３に示すように、尾根谷度は樹高データをもとに計算した地上開度と地下開度から算出される。このため、樹冠が隣木に接していない樹木は樹冠の端でDHM が急激に変化し、尾根谷度が大きくなる。これにより、樹冠端部が樹木頂点より大きな尾根谷度になる。
【００２１】
図４は尾根谷度画像と本実施の形態２の樹冠形状指数との比較を説明する説明図であり、図４（ａ）は尾根谷度が樹冠端部において極大化（画像が明るい領域）している例を示す画像である。
【００２２】
この極大化により、林縁部では、樹冠端部が頂点として認識されてしまったり、樹木頂点が樹冠内に複数生じたりする可能性がある。この問題を解決するには、樹冠端部の極大化を防ぎ、樹木の頂点の尾根谷度を強調して隣木との区別を明確にする方法が好ましい。
【００２３】
図５は実施の形態の樹木頂点認識装置の概略構成図である。図６は本実施の形態を説明する概略フローチャートである。実施の形態の樹木頂点認識装置の概略構成図を説明する前に、図６の概略フローを説明する。図６に示すように、航空レーザ計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ２０）。また、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ２１）。
【００２４】
そして、DSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）を樹木高グリッドデータＤｄｉとして求めて（Ｓ２２）、樹冠形状の再現を行う（Ｓ２３）。
【００２５】
樹冠形状の再現は、この樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）を用いて例えば５ｍ範囲（図７を参照）における各グリッドデータと自身のデータとのアスペクト比を求め、このアスペクト比を一定条件下で新たなアスペクト比に変更し、変更したアスペクト比を用いて地上開度φ１´及び地下開φ２´度を算出し（Ｓ２３ａ）、図４（ｂ）に示す樹冠形状指数画像のデータを算出する（Ｓ２３ｂ）。この樹冠形状指数画像については後述する。
【００２６】
そして、この樹冠形状指数Gｅiから樹木の頂点を抽出する（Ｓ２４）。樹木頂点の抽出は、樹冠形状指数Gｅiに局所最大値フィルタ（3×3範囲）を適用して（Ｓ２４ａ）、「樹木頂点の候補点」を抽出する（Ｓ２４ｂ）。
【００２７】
また、本実施の形態では、後述する動的閾値により樹冠部を抽出して（Ｓ２４ｃ）、樹冠部画像を得る（Ｓ２４ｄ）。
【００２８】
そして、頂点候補が存在する樹冠部領域毎に樹木頂点候補の中から最高の樹木高を有する点を求め、これを樹木頂点とし（Ｓ２４ｅ）、抽出する（Ｓ２４ｆ）。これによって得られた画像を図８に示す。
【００２９】
次に、図５の構成を説明する。
【００３０】
（図５の説明）
図５においては、表示処理部１９は、表示部１３（画面）に各データベース、メモリ等の画像を表示させたり、画面の画像データを読み込む等の機能を有しているが、各データベース、メモリ等とのやり取りを示す線については記載を省力している。
【００３１】
また、図５においては、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を記憶したデータベース１０と、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model）を記憶したデータベース１１と、特定地域のデジタル写真（オルソフォト画像）を記憶したデータベース１２と、樹木高算出部１４と、樹木高データが記憶されるファイル（メモリ）１５と等を備えている。
【００３２】
また、検索範囲決定部１７と、表示処理部１９と、局所最大フィルタ部２０と、樹木頂点候補抽出部２１と、樹木頂点候補データが記憶されるファイル２２と、樹冠部抽出部２３と、樹冠部のデータが記憶されるファイル２４と、第１の樹木頂点決定部２５と、樹木頂点データが記憶されるファイル２９と出力部３０等を備えている。
【００３３】
さらに、本実施の形態の樹木頂点認識装置は、第１の樹冠形状指数算出部３２と、閾値算出部３５と、局所領域設定部３６と等を備えている。
【００３４】
データベース１０は、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model：本実施の形態では０．５ｍ間隔の格子状のデータ）を記憶している。このグリッドデータは、地盤のＸ、Ｙ、Ｚ（地盤高Ｚａ）とを対応させて記憶させられている。つまり、レーザデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて地面を復元したものである。
【００３５】
データベース１１は、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model：本実施の形態では０．５ｍ間隔）を記憶している。このグリッドデータ（本実施の形態では０．５ｍ間隔の格子状のデータ）は、樹木表面のＸ、Ｙ、Ｚ（Ｚｂ：樹木表面高）とを対応させて記憶している。つまり、レーザデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて表層を復元したものである。
【００３６】
すなわち、これらのグリッドデータは、格子間隔０．５ｍのＸ−Ｙ平面直角座標系に、緯度経度に対応するＸ値、Ｙ値とレーザ計測によって得られたＺ値を割り当てている。
【００３７】
樹木高算出部１４は、データベース１０のＤＳＭとデータベース１１のＤＥＭとの差を求め、これをファイル１５（メモリ）のＸ，Ｙ座標に割付ける（以下樹木高グリッドデータＤｄｉという）。
【００３８】
検索範囲決定部１７は、表示部１３の画面上においてオペレータによって決められた検索範囲Ｍｉを樹冠形状指数算出部３２に設定する。
【００３９】
第１の樹冠形状指数算出部３２は、デジタル写真を基に決定された検索範囲Ｍｉに対応する樹木高グリッドデータＤｄｉをファイル１５から読み込んで、これらのグリッドデータＤｄｉにおけるアスペクト比（樹木高Ｚｃｉ/距離ｄｉ）の最大、最小値を基にして樹冠形状を明瞭にする所定の値に変更する（第１のアスペクト比更新処理３１ａ）。
【００４０】
そして、この変更された新たな最大値、新たな最小値を用いて新たな地上開度φ１´（以下第１の地上開度φ１´という）、地下開度φ２´（以下第１の地下開度φ２´という）を求め、これらの開度から算出される第１−１の樹冠形状指数φ３´を、例えば８方向で平均化して、平均化した樹冠形状指数ｅｉ（以下第１−２の樹冠形状指数ｅｉという）を得る（第１の樹冠形状指数計算処理３１ｂ）。
【００４１】
次に、局所領域Ｊｉ(例えば3×3pixel)で平滑化し、これを検索範囲Ｍｉ（例えば５ｍ）における第１−３の樹冠形状指数Ｇｅｉとしてファイル（メモリ）３７に保存する（平滑化処理３１ｃ）。
【００４２】
前述の第１−３の樹冠形状指数Ｇｅｉは、設定された樹木高グリッドデータの番号と、Ｘ、Ｙ座標と、樹木高の差Ｚと、第１−１の樹冠形状指数φ３´と、第１−２の樹冠形状指数ｅｉと、この第１−１の樹冠形状指数φ３´を得たときの地上開度φ１´と、地下開度φ２´と、検索方向毎の検索範囲以内の樹木高グリッドデータ（Ｘ、Ｙ、Ｚ）と、これらの樹木高グリッドデータまでの距離ｄｉと、設定された樹木高グリッドデータに対して高さの差Ｚ等にリンク付けされている。これらを総称して第１の樹冠形状指数情報Ｇｉという。
【００４３】
閾値算出部３５は、表示部１３に表示された画像（グリッドデータ又はオルソフォト画像又はオルソフォト画像とグリッドデータとを重ね表示）内で、予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第１−３の樹冠形状指数Ｇｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【００４４】
局所最大値フィルタ２０は、ファイル３７の第１の樹冠形状指数情報Ｇｉの各樹木高グリッドデータ（Ｘ、Ｙ、Ｚ）に対応する第１−３の樹冠形状指数Ｇｅｉに局所領域Ｆｉ（例えば３×３範囲）のフィルタを順次適用し、この結果のデータＤｆｉをファイル２７に保存する。
【００４５】
樹木頂点候補抽出部２１は、ファイル２７のデータＤｆｉとファイル３７の第１−３の樹冠形状指数Ｇｅｉとを比較し、同じ値のものを頂点候補Ｄｇｉとしてファイル２２に保存する。
【００４６】
樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の、第１−３の樹冠形状指数Ｇｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル２４に保存する。
【００４７】
第１の樹木頂点決定部２５は、ファイル２４の樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ、一致したものを仮の頂点と決定し、これをファイル２８に保存する。そして、同一樹冠部に複数の仮の頂点が存在するときは、この樹冠部で最も樹木高が高いものを単木の頂点として決定し、これをファイル２９に保存する（Ｐｉ）。
【００４８】
表示処理部１９は決定した頂点を例えば○印（△印、×印、番号等の識別）にして画面に表示する。このとき、写真画像と重ねて表示してもよい。
【００４９】
出力部３０は、ファイル２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。このとき、樹木高も表示してもよい。
【００５０】
上記のように構成された実施の形態１の樹木頂点認識装置について以下にフローチャートを用いて説明する。
【００５１】
図９は樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【００５２】
本実施の形態では、初めに初期設定を行う（Ｓ３０）。この初期設定は、例えば画面のオルソフォト画像において所定のエリアの設定、検索範囲Ｍｉ、検索方向ｋｉ等の設定である。
【００５３】
前述の検索方向ｋｉは、以下のようにして求めている。図9のS32において、任意のDdipを決定したときに、隣接するDdiを検索し、データが存在する方向の数を検索方向kiとする。例えば、任意のDdipが画像の右上角であった場合、隣接データ数は3となり、検索方向は3となる。また、任意のDdipが画像の中央であった場合、隣接データ数は8となり、検索方向は8となる。
【００５４】
本実施の形態では検索範囲Ｍｉはオルソ画像を確認し、対象林分の代表的な樹木間隔未満で設定する。
【００５５】
そして、樹冠形状指数算出部３２は、検索範囲の樹木高グリッドデータＤｄｉ（DSMとDEMデータの差分）をファイル１５から全てリードする（Ｓ３１）。
【００５６】
次に、これらの樹木高グリッドデータＤｄｉの中の１つの樹木高グリッドデータＤｄｉｐの番号を設定すると共に、検索方向内の１つの方向ｋｉｐを設定する（Ｓ３２、Ｓ３３）。本実施の形態では８方向として説明する。
【００５７】
次に、設定された番号の樹木高グリッドデータＤｄｉｐから方向ｋｉｐで検索範囲Ｍｉ内の樹木高グリッドデータＤｄｉｐ（以下グリッドデータＤｐｐｉという）を抽出する（Ｓ３４）。
【００５８】
前述の地上・地下開度を算出する検索範囲Ｍｉは対象林分に応じて変化させる必要がある。検索範囲が樹木間隔より長過ぎると隣木の頂点に対して最大傾斜角をなすため、自身が頂点として検出できない可能性がある。また、樹木間隔に対して短過ぎると樹冠内の山谷形状を検出し、一つの樹冠に対して複数の樹木頂点が抽出される場合がある。以上から検索範囲は対象林分の代表的な樹木間隔未満で設定するのが良い。
【００５９】
そして、設定された番号の樹木高グリッドデータＤｄｉｐからグリッドデータＤｐｐｉまでの間の距離ｄｉを求める（Ｓ３５）。本実施の形態では、例えば5mとして説明する。
【００６０】
次に、設定された番号の樹木高グリッドデータＤｄｉｐとグリッドデータＤｐｐｉとの垂直方向の樹木高差ΔＺｃｉ（差であることを強調するためにΔを用いている）を求める（Ｓ３６）。
【００６１】
次に、この樹木高差ΔＺｃｉを距離ｄｉで割った値（アスペクト比Ａｄｉ）を求め（Ｓ３７）、これを検索方向ｋｉｐと樹木高グリッドデータＤｄｉｐと、この番号とに対応（総称してグリッドデータ毎のアスペクト比データＶｉという）させてファイル３７に記憶する（Ｓ３８）。
【００６２】
次に、図１０に示すように、検索方向ｋｉｐが奇数か偶数かどうかを判定する（Ｓ３９ａ）。
【００６３】
ステップＳ３９ａにおいて、検索方向ｋｉｐが偶数と判定したときは、グリッドデータＤｐｐｉが７個かどうかを判定する（Ｓ３９ｂ）。ステップＳ３９ｂにおいて、グリッドデータＤｐｐｉが７個と判定したときは、処理を後述するステップＳ４１に移す。
【００６４】
また、ステップＳ３９ｂにおいて、７個ではないと判定したときは、グリッドデータＤｐｐｉを次のグリッドデータに更新して処理を図９のステップＳ３４に戻す（Ｓ３９ｃ）。
【００６５】
また、ステップＳ３９ａにおいて、検索方向ｋｉｐが奇数と判定したときは、グリッドデータＤｐｐｉが１０個かどうかを判定する（Ｓ３９ｄ）。
【００６６】
ステップＳ３９ｄにおいて、１０個と判定したときは、処理を後述するステップＳ４１に移す。
【００６７】
前述のステップＳ３９ａ、３９ｂ、３９dの判定は、例えば０．５ｍメッシュにおいて検索範囲Ｍｉが５ｍであれば、上下左右方向の場合（奇数の場合）は１０個入るので１０個かどうかを判定している。また、斜め方向の場合（偶数の場合）は７個入るので７個かどうかを判定している。
【００６８】
図１１においては、５ｍではなく１ｍの検索範囲を示して、検索範囲Ｍｉ内の検索個数を説明している。
【００６９】
また、ステップＳ４１においては、検索方向ｋｉｐが８方向になったかどうかを判定する（Ｓ４１）。
【００７０】
ステップＳ４１において、検索方向ｋｉｐが８方向全てになっていないと判定したときは検索方向を次の検索方向に更新して（Ｓ４２）、処理をステップＳ３３に戻す。
【００７１】
また、ステップＳ４１において検索方向ｋｉｐが８方向全てと判定したときは、エリア内の樹木高グリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）が最後かどうかを判定し（Ｓ４３）、最後でないときは、樹木高グリッドデータＤｄｉの番号を更新して（Ｓ４４）、処理を図９のステップＳ３２に戻す。
【００７２】
また、ステップＳ４３において、樹木高グリッドデータＤｄｉが最後と判定したときは、アスペクト比更新処理を起動する（Ｓ４５）。
【００７３】
このアスペクト比更新処理の方法の前にアスペクト比と第１の地上・地下開度（φ１´、φ２´）との関係を説明する。
【００７４】
地上・地下開度の算出には任意の点から比較する点の距離と標高差が必要となる。算出式は次式で表される。
【００７５】
Max_angle = tan-1 (最大ΔZ/d)・・・・・・式（A）
Min_angle = tan-1 (最小ΔZ/d)・・・・・・式（B）
ここで、ΔZは標高差、dは2点間の距離を表す。
【００７６】
第１の地上開度φ1´= 90−Max_angle・・・・・・式（C）
第２の地下開度φ2´= 90+ Min_angle・・・・・式（D）
方向毎に求めたアスペクト比（上記例では１０個）の中から最大及び最小のアスペクト比を抽出し、式Aと式Bから最大傾斜角と最小傾斜角を求め、式Cと式Dに当てはめて地上・地下開度（φ１´、φ２´）を算出する。但し、最大傾斜角と最小傾斜角の計算は後述するアスペクト比の更新の後となる。
【００７７】
第１の地上開度φ１´及び第２の地下開度φ２´を求めるためのアスペクト比は例えば次のように更新する（図１２参照）。
【００７８】
（1）.「ΔZ/d」が−1.7320（−６０度）以下の場合は-1.7320（-60 度）とする。
【００７９】
（樹冠端部の極大化防止）
（2）.「ΔZ/d」が-1.7320（-60 度）から0（0 度）の場合は-1.7320（-60 度）とする。
【００８０】
（樹冠山地形の強調）
（3）.「ΔZ/d」が0（0 度）から1.7320（60 度）の場合は1.7320（60 度）とする。
【００８１】
（樹冠谷地形の強調）
（4）.「ΔZ/d」が1.7320（60 度）以上の場合は「ΔZ/d」の値を据え置く。
【００８２】
（5）.「ΔZ/d」が隣接する8 方位全てがマイナスである場合「ΔZ/d」を1000倍する。
【００８３】
（樹木頂点の強調）
という処理である。
【００８４】
前述の（1）から（4）で60度を基準としてアスペクト比を更新しているが、これは現地調査でスギの樹冠形状が60度程度であったため、現地の樹冠形状に適した値として採用した。
【００８５】
すなわち、ここでは６０度と設定しているが、想定する樹冠形状によりオペレータが任意に角度を設定する。
【００８６】
更に、（５）では樹木頂点を強調するために「ΔZ/d」を1000 倍している。これにより最小角はほぼ-90度となり、最大限の樹木頂点強調が出来る。具体的には、（1）、（2）から「ΔZ/d」の変更後の値は-1.7320 となり、これを1000倍して角度を算出すると-89.9度となる。
【００８７】
このように算出された第１の地上・地下開度（φ１´、φ２´）を用いて尾根谷度の算出式（φ3 =(φ1 −φ2)/2）
と同様にして第１−１の樹冠形状指数φ3´を求める。この第１−１の樹冠形状指数φ3´を8 方位で計算し、平均化したものを第１−２の樹冠形状指数ｅｉとする。
【００８８】
以下に第１−２の樹冠形状指数ｅｉを求めるまでの処理を具体的に説明する。初めにアスペクト比更新処理について説明する。
【００８９】
アスペクト比更新処理について図１３のフローチャートを用いて説明を補充する。図１３に示すように、アスペクト比更新処理においては、ファイル３７の樹木高グリッドデータＤｄｉを設定し、検索方向ｋｉｐを設定する（Ｓ５０、Ｓ５１）。
【００９０】
次に、設定された番号の樹木高グリッドデータＤｄｉの検索方向ｋｉｐ内でのアスペクト比データViの中から値が最大なアスペクト比（ＭＡＸＡｄｉ）、値が最小のアスペクト比（ＭＩＮＡｄｉ）を抽出（読込む）する（Ｓ５２）。
【００９１】
例えば、0.5mグリッドデータを用いて、検索範囲Ｍｉが1mであるとき、ある検索方向において、アスペクト比が-2m/0.5mと-2m/1.0mであったとすると、これらの比は-4.0と-2.0になり、最大値が-2.0、最小値が-4.0となる。
【００９２】
次に、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉが上記（３）の０〜１．７３２０の範囲（角度では０度〜６０度）かどうかを判定する（Ｓ５３）。
【００９３】
ステップＳ５３において、０〜１．７３２０の範囲と判定したときは、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉを１．７３２０（ＭＡＸＡｄｉ、ＭＩＮＡｄｉ：Ａｄｂｉ）として処理をステップＳ５９に移す（Ｓ５４）。これにより、頂点から下の場合は、頂点から下の谷地形であることを強調する。
【００９４】
また、ステップＳ５３において、０〜１．７３２０の範囲（角度では０度〜６０度）ではないと判定したときは、読み込んだアスペクト比が上記（４）の１．７３２０以上かどうかを判定する（Ｓ５５）。
【００９５】
ステップＳ５５において、アスペクト比が１．７３２０以上と判定したときはそのアスペクト比（Ａｄｉ）をそのままの値（Ａｄａｉ）にして処理をステップＳ５９に移す（Ｓ５６）。
【００９６】
また、ステップＳ５５において、読み込んだアスペクト比が１．７３２０以上ではないと判定したときは、読み込んだアスペクト比が上記の（１）の−１．７３２０以下かどうかを判定する（Ｓ５７）。
【００９７】
ステップＳ５７において、−１．７３２０以下であると判定したときは、ステップＳ５９-１に処理を移す。
【００９８】
ステップＳ５９−１においては、抽出したアスペクト比を−１．７３２０（Ａｄｃｉ）にする。
【００９９】
これにより、樹冠端部と地盤との位置関係から生じるアスペクト比の極小化（樹冠形状指数では極大化）を防止する。
【０１００】
また、ステップＳ５７において、−１．７３２０以下ではないと判定したときは、抽出したアスペクト比（０〜−１．７３２０（２）(０を除く)）を−１．７３２０（Ａｄｃｉ）にする（Ｓ５９−２）。
【０１０１】
これにより、樹冠端部から上の場合は、樹冠端部から上の山地形であることを強調する。
【０１０２】
また、前述のステップＳ５４、ステップＳ５６、ステップＳ５９−２、ステップＳ５９−１のアスペクト比をファイル３７に記憶する（Ｓ５９）
次に、検索方向ｋｉｐが８かどうかを判定し（Ｓ６１）、８方向になっていないときは検索方向を更新して処理をステップＳ５１に戻す（Ｓ６２−１）。
【０１０３】
ステップＳ６２において、隣接する８方向の樹木高グリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスと判定したときは、８方向の更新後のアスペクト比を１０００倍（Ａｄｄｉ）する（Ｓ６３）。すなわち、頂点を強調する。
【０１０４】
次に、検索エリア内で樹木高グリッドデータＤｄｉが最後かどうかを判定し（Ｓ６４）、最後でないときは次の番号の樹木高グリッドデータＤｄｉに更新して処理をステップＳ５０に戻す（Ｓ６５）。
【０１０５】
すなわち、ファイル３７には、ステップＳ５９の処理によって、図１４に示すように、設定された番号の樹木高グリッドデータＤｄｉ（図１４においてはｉ＝１）と検索方向ｋｐｉと抽出したＭＡＸＡｄｉ、ＭＩＮＡｄｉと更新したアスペクト比とが対応させられて保存されている。
【０１０６】
更新したアスペクト比は、アスペクト比Ａｄｉのままの値（Ａｄａｉ）又は１．７３２０（Ａｄｂｉ）又は−１．７３２０（Ａｄｃｉ）若しくは上記3種類のいずれかの値を１０００倍したアスペクト比（Ａｄｄｉ）のいずれかである。
【０１０７】
これらを総称して本実施の形態では更新されたアスペクト比を更新アスペクト比と称する。
【０１０８】
次に、第１の樹冠形状指数算出部３２の樹冠形状指数の算出を図１６のフローチャートを用いて説明する。検索エリア内の各樹木高グリッドデータからの検索範囲において全てのアスペクト比の更新が行われると、図１６に示すように樹木高グリッドデータＤｄｉの番号を設定し（Ｓ７０）、検索方向ｋｉｐを設定する（Ｓ７１）。
【０１０９】
次に、設定されたＤｄｉ、ｋｉｐに対応する更新アスペクト比を引当、この更新アスペクト比の中から最大アスペクト比ＭＡＸＡｄｉ（アスペクト比Ａｄｉのまま又は１．７３２０又は−１．７３２０又は前記Ａｄｉを１０００倍したアスペクト比又は1732又は-1732のいずれか）と、最小アスペクト比ＭＩＮＡｄｉ（アスペクト比Ａｄｉのまま又は１．７３２０又は−１．７３２０又は前記Ａｄｉを１０００倍したアスペクト比又は1732又は-1732のいずれか）読み込む（Ｓ７２）。
【０１１０】
そして、第１の地上開度φ１´と第１の地下開度φ２´を算出する。（Ｓ７３）。
【０１１１】
第１の地上開度φ１´＝９０度−tan-1(ＭＡＸＡｄｉ)
第１の地下開度φ２´＝９０度＋tan-1(ＭＩＮＡｄｉ)
次に、この第１の地上開度φ１´と第１の地下開度φ２´とを用いて第１−１の樹冠形状指数φ３´を求める（Ｓ７４）。
【０１１２】
第１−１の樹冠形状指数φ３´＝（φ１´−φ２´）／２
次に、８方位全てに対して樹冠形状指数を求めたかどうかを判断する（Ｓ７５）。ステップＳ７５において８方位全てに対して第１−１の樹冠形状指数φ３´を求めていないと判定したときは、検索方向を更新して処理をステップＳ７１に戻す（Ｓ７６）。
【０１１３】
次に、ステップＳ７５で８方位全てに対して第１−１の樹冠形状指数φ３´を求めたと判定したときは、これを平均化（８方位）する（Ｓ７７）。
【０１１４】
そして、この平均化された第１−２の樹冠形状指数ｅｉを、樹木高グリッドデータＤｄｉを中心とした検索範囲での第１−２の樹冠形状指数ｅｉとして（Ｓ７８）、樹冠グリッドデータＤｄｉにリンク付けしてファイル３７に記憶する（Ｓ７９）。
【０１１５】
次に、検索エリア内の樹冠グリッドデータＤｄｉの番号が最後かどうかを判定し、最後でないときは樹冠グリッドデータＤｄｉの番号を更新して処理をステップＳ７０に戻す。
【０１１６】
最後に第１−２の樹冠形状指数ｅｉを局所領域Ｊｉ（例えば3×3pixel）で平滑化したものを第１−３の樹冠形状指数Ｇｅｉとして検索範囲に対応させてファイル３７に記憶する。
【０１１７】
すなわち、ファイル３７には図１４に示すように、設定された樹木高グリッドデータＤｄｉと、検索方向ｋｉｐと、更新前のＭＡＸ、ＭＩＮのアスペクト比と、更新後のＭＡＸ、ＭＩＮの更新アスペクト比とが検索方向毎に、リンクつけされて記憶されている。
【０１１８】
また、ファイル３７には図１５に示すように、設定された樹木高グリッドデータＤｄｉの番号と、そのＸ、Ｙ座標と、樹木高の差Ｚと、次の番号の樹木高グリッドデータと、このＸ、Ｙ座標と、距離ｄｉと、第１−１の樹冠形状指数φ３´と、この第１−１の樹冠形状指数φ３´を得たときの地上開度φ１´と、地下開度φ２´とが検索方向にリンクつけされ、また、８方向全てで平均化した第１−２の樹冠形状指数ｅｉ等が対応させて記憶されている。図１５には、平滑後の第１−３の樹冠形状指数Ｇｅｉを対応させて例示している。
【０１１９】
次に、頂点候補の抽出、閾値の算出方法、樹冠部の抽出処理の概略を以下に説明する。
【０１２０】
頂点候補の抽出
ここでは樹冠形状指数から樹木頂点候補を抽出する。抽出方法は次の通りである。樹冠形状指数に局所最大値フィルタを適用し、局所領域Ｆｉ内の最大値を中央画素値に置き換える。この際、局所領域Ｆｉを例えば3×3に設定した。このフィルタ処理後の指標と処理前の指標を比較し、処理前後で第１−３の樹冠形状指数Ｇｅｉが同値であったときに、その地点を「樹木頂点候補」として抽出する。
【０１２１】
局所最大値フィルタ２０は、第１−３の樹冠形状指数Ｇｅｉに対して、局所領域Fiの局所最大値フィルタを適用して、これをファイル２７に保存（Ｄｆｉ）する。
【０１２２】
次に、樹木頂点候補抽出部２１がファイル２７のデータＤｆｉを頂点候補（以下Ｄｇｉ）として抽出し、ファイル２２に保存する。具体的にはファイル３７のデータの値とファイル２７の同じ座標値のデータの値が同値の場合を「１」と、他を「０」としてファイル２２に定義する（図２０参照）。
【０１２３】
ここで、局所最大フィルタ部２０及び樹木頂点候補抽出部２１の処理について説明を補充する。
【０１２４】
局所最大フィルタ部２０（Local max filter２０ともいう）は、ファイル３７の第１の樹冠形状情報Ｇｉの第１−３の樹冠形状指数Ｇｅｉ（Ｇｅ１、Ｇｅ２・・）に局所領域Fi内で最大となる値を局所領域Fiの中央画素値に置き換える処理を順次行う。これらをＤｆｉとしてファイル２７に保存する。
【０１２５】
そして、樹木頂点候補算出部２１がデータＤｆｉとファイル３７のＧｅｉとを比較し、値が同じものを頂点候補Ｄｇｉとしてファイル２２に保存している。
【０１２６】
次に樹冠部を抽出する閾値の算出方法及び樹冠部の抽出処理を説明する。
【０１２７】
閾値以上の領域を樹冠部として抽出するが、これまで閾値は対象とする林分のオルソ画像で単木を目視確認して、同様に樹木が分割される値を閾値としていた。この手法で設定される閾値は一定値であることから樹冠部であっても閾値より小さい場合は樹冠部として抽出できなかった。この問題を解決するため、新たに画像に応じた閾値を局所領域Ｈｉ（3×3pixel や5×5pixel など）による算出で動的に決定する。
【０１２８】
初めに一つの樹冠領域を想定する。更に、この一つの樹冠領域に対して樹冠部として抽
出する領域を想定する（図１７）。ここでは尾根形状による隣木との同一化を防ぐために、樹冠領域に対して半分程度を樹木頂点が一つ含まれる樹冠部として抽出する。樹冠領域を3×3pixel から9×9pixel を想定すると樹冠部は1×1 から4×4 となり、樹冠領域に対する樹冠部の割合は平均して16％となる（表１）。
【表１】
【０１２９】
なお、0.5m グリッドデータ使用していることから、樹冠領域は1.5m×1.5m から4.5m×4.5m となり、現地で確認した樹冠領域と同等と考える。この樹冠領域を局所領域として、樹冠部の割合が16％程度になるように、局所領域の平均値（μ）に標準偏差（σ）を足した値を閾値とした。設定根拠は次の通りである。正規分布しているヒストグラムにおいて平均値の±σの範囲は全pixel の68％を占める。更に最小値から平均値の＋σの範囲を考えると84％に該当し、μ＋σ以上の範囲が16%を占める（図１８）。
【０１３０】
具体的には局所領域Ｈｉ内の平均値と標準偏差を求め、平均値と標準偏差の和を中央画素値に置き換える。第１−３の樹冠形状指数Ｇｅｉと動的閾値Ｈｇｉを比較して樹冠形状指数が大きければ、樹冠部として抽出する。
【０１３１】
次に、図１９のフローチャートを用いて具体的に樹冠部の抽出処理を説明する。初めに局所領域Ｈｉを設定する（Ｓ８０）。
【０１３２】
本実施の形態では所定エリアの樹木高グリッドデータＤｉが画像化されてオルソフォト画像に重ね表示されているとする。画面に表示された樹木高グリッドデータＤｄｉを本実施の形態では画像化されたグリッドデータＲｉと称する。
【０１３３】
この状態で、オペレータは例えば３ピクセル×３ピクセルの局所領域Ｈｉを入力する。樹冠部決定処理は、局所領域、検索エリアをメモリ（図示せず）に設定する初期設定を行う（Ｓ８０）。
【０１３４】
次に、閾値算出部が検索エリアの画像化されたグリッドデータＲｄｉを全てリードする（Ｓ８１）。
【０１３５】
次に、これらのグリッドデータＲｄｉの中からグリッドデータＲｄｎを設定し（Ｓ８２）、設定した局所領域Ｈｉに基づき、グリッドデータＲｄｎ（Ｘｄｎ、Ｙｄｎ、樹木高dn）を中心に、このグリッドデータの座標値を有する第１の樹冠形状情報Ｇｉ（３×３）をファイル３７からリードする（Ｓ８３）。
【０１３６】
このリードされた３×３の第１の樹冠形状情報Ｇｉは、第１−３の樹冠形状指数ＧｅｉをＺｉとし、グリッドデータＤｄｉのＸ軸をＤｘｉ、Ｙ軸をＤｙｉと記載する。
【０１３７】
例えば、
（Ｄｘ１、Ｄｙ１、Ｚ１）、（Ｄｘ２、Ｄｙ２、Ｚ２）、（Ｄｘ３、Ｄｙ３、Ｚ３）
（Ｄｘ４、Ｄｙ４、Ｚ４）、（Ｄｘ５、Ｄｙ５、Ｚ５）、（Ｄｘ６、Ｄｙ６、Ｚ６）
（Ｄｘ７、Ｄｙ７、Ｚ７）、（Ｄｘ８、Ｄｙ８、Ｚ８）、（Ｄｘ９、Ｄｙ９、Ｚ９）
Ｚｉ：平均化、平滑された第１−３の樹冠形状指数Ｇｅｉ
となる。
【０１３８】
そして、閾値算出部３５がこれらの第１−３の樹冠形状指数Ｇｅｉの平均値（μ）と標準偏差（σ）とを算出し（Ｓ８４）、これらの和を求め、これを設定した局所領域Ｈｉの閾値（Ｈｇｉ）として求める（Ｓ８５）。
【０１３９】
次に、閾値算出部３５は、設定したＲｄｎに対応した閾値を樹冠部抽出部２３に出力する（Ｓ８６）。樹冠部抽出部２３は第１−３の樹冠形状指数Ｇｅｉとこの閾値とを比較する（Ｓ８７）。
【０１４０】
そして、閾値（Ｈｇｉ）に対して第１−３の樹冠形状指数Ｇｅｉが大きいかどうかを判定する（Ｓ８８）。ステップＳ８８において、第１−３の樹冠形状指数Ｇｅｉの方が大きいと判定したときは、樹冠部とする（Ｓ８９）。
【０１４１】
次に、グリッドデータＲｄｉは最後かどうかを判定し（Ｓ９１）、最後でないときはグリッドデータＲｄｉを次のグリッドデータに更新し、処理をステップＳ８２に戻す（Ｓ９２）。
【０１４２】
また、ステップＳ８８において閾値（Ｈｇｉ）の方が大きいと判定したときは、樹冠部ではないと判定する（Ｓ９０）。
【０１４３】
また、樹冠部抽出部２３が閾値Ｈｇによって樹冠部のデータＤｈｉを抽出し、これをファイル２４に保存している（図２１参照）。
【０１４４】
図２１（ａ）のデータＤｈｉは、閾値を超えるＧｅｉを抽出したものであり、例えば番号１のエリアの第１−３の樹冠形状指数Ｇｅｉは、（Ｘｈ５、Ｙｈ５）を「１」、（Ｘｈ６、Ｙｈ６）を「１」・・（Ｘｈ１０、Ｙｈ１０）を「１」として示している。
【０１４５】
図２１の（ｂ）は図２１（ａ）を具体的に示すものであり、閾値を超えた画素は「１」とする。「１」が連続している範囲が樹冠部である。
【０１４６】
そして、ファイル２４の中の頂点候補ＤｈｉのＸｈｉ、Ｙｈｉをポリゴン化（座標Ｘｈｉ、Ｙｈｉと色情報）して表示処理部１９によってラベルリング（纏まったエリアを（樹冠部）にラベルを振る）させる（図２６参照）
次に、樹木頂点決定部２５は、ファイル２４のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定する。そして、樹冠部のエリア内に複数の仮の頂点が存在するときは、この樹冠部のエリアで最も樹木高が高いものを単木の頂点として決定し、これをファイル２９に保存する（Ｐｉ）。表示処理部１９は決定した頂点を画面に表示する（図２３参照）
前述の頂点の決定について説明を補充する。例えば、図２２の（ａ）に示すように、ファイル２４に第１列目が「００１１００１１１００」、第２列目が「０１１１１０１１１１０」、第３列目が「００１１０００１０００」、第４列目が「０００００１０００１１」という具合に頂点上部域が得られた場合に、図２２（ｂ）に示すように、「１」が隣接するまとまり（「０」が区切り）を１つの樹冠部としてまとまり毎に通し番号を割り振っていく。
【０１４７】
そして、例えば、同一高の頂点が樹冠内に複数の頂点が存在することになると、一度、ファイル２８のデータＤＱｉを表示させて、オペレータによって指定させ、指定されたものを頂点Ｐｉと決定してファイル２９に保存する。
【０１４８】
このような処理によって第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出結果は、図８、図２３に示すように頂点が求められている。
【０１４９】
また、表２に尾根谷度と第１−３の樹冠形状指数Ｇｅｉの樹木頂点抽出結果の比較を示す。
【表２】
【０１５０】
表２から、抽出精度は同程度であるが、抽出率及び正解率が向上していることが分かる。以上から第１−３の樹冠形状指数Ｇｅｉの有意性が確認できた。
【０１５１】
また、図２４の尾根谷度画像では林道に接した樹木は極大化して画像が明るくなっているが図２５の樹冠形状指数画像では同箇所が暗くなり、極大化を解消していることがわかる。
【０１５２】
次に、尾根谷度及び第１−３の樹冠形状指数による樹木頂点の抽出数、判読による樹木頂点数と現地で樹冠上層部に到達していると確認した樹木本数を表３にまとめる。
【表３】
【０１５３】
表から、尾根谷度による樹木頂点抽出はP8-2 を除いて3〜4 割の程度の抽出率であったのに対して、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出は7〜9 割の抽出率であった。また、空中写真からの情報では密度の高い森林において樹木の抽出に限界があることがわかっている。このことから、同様に上空からの情報であるレーザデータを用いた樹木頂点抽出では判読による樹木頂点数を検証値として考えることもできる。
【０１５４】
そこで、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出と判読による樹木頂点の割合を算出した。その結果、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出数は判読の樹木頂点数の8 割強〜10 割強であった。以上から、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出は尾根谷度による抽出よりも倍近い抽出率が得られ、なおかつ、目視判読による樹木頂点抽出とほぼ同レベルで行えることが検証できた。
【０１５５】
また、本実施の形態は、図２７に示すようなコンピュータシステムVPSによって実現するのが好ましい。ＶＰＳ１は、ワークステーション、プロセッサ、マイクロコンピュータ、ロジック、レジスタ等の適宜な組み合わせからなる中央情報処理装置（CPU）と、この中央情報処理装置に必要な制御・操作情報を入力するキーボード（KB）、マウス、対話型ソフトスイッチ、外部通信チャネル等を含む情報入力部と、中央情報処理装置から出力された情報を広義な意味で表示・伝送するディスプレイ、プリンタ、外部通信チャネル等を含む情報出力部と、中央情報処理装置に読み込まれるオペレーティングシステム、アプリケーションプログラム等の情報が格納されたロム（ＲＯＭ）と、中央情報処理装置で随時処理すべき情報及び中央情報処理装置から随時書き込まれる情報を格納するラム（ＲＡＭ）等とを備える。ＲＯＭ、ＲＡＭを適宜統合、細分化することは差し支えない。
【０１５６】
前述のアプリケーションプログラムは、本実施の形態の、樹木高算出部１４と、検索範囲決定部１７と、局所領域設定部３６と、第１の樹冠形状指数算出部３２と、閾値算出部３５と、局所最大フィルタ２０（Ｌｏｃａｌ max ｆｉｌｔｅｒともいう）と、樹木頂点候補抽出部２１と、樹冠部抽出部２３と、第１の樹木頂点決定部２５と、出力部３０等である。
【０１５７】
なお、上記実施の形態では、画面において決定した検索エリアにおける樹木の各樹木頂点を求めるとしたが、検索エリアではなく、樹木画像上で軌跡Ａｉを描き、この軌跡Ａｉ上における樹木の頂点を抽出してもよい。
【０１５８】
つまり、軌跡Ａｉの樹木高グリッドデータを抽出し、この樹木高グリッドデータに対応する第１−３の樹木形状指数Ｇｅｉを読み出し、この樹木形状指数Ｇｅｉによって頂点を決定する。
【０１５９】
＜実施の形態２＞
スギ、ヒノキ人工林は通常3000 本/ha 植栽され、林齢17 年程度で第1 回間伐が行われる。伐採率は3 割程度であるが、自然枯死を併せると間伐後は1950 本/ha 程度となる。
【０１６０】
その後、6〜8 年毎に第2 回、第3 回の間伐（間伐率は共に約3 割）を行い、950 本/ha 程度の立木密度にして収穫を行う。現在日本のスギ、ヒノキの植林は少なく、そのほとんどは林齢20 年以上であり、第1 回の間伐は実施されていることと想定される。
【０１６１】
このことから、第2 回の間伐が終了しているか判断するためにも2000 本/ha 程度の立木密度の森林を把握することが求められる。
【０１６２】
本実施の形態２は2000 本/ha 程度の立木密度の密な森林（例えばスギ、ヒノキ等）においても樹木頂点の抽出精度を向上できる。
【０１６３】
本実施の形態２では、抽出精度の向上のために、「地上開度、地下開度の算出の見直し」と「局所加重平均」と、「樹木頂点が隣接する場合に所定の樹木頂点を除去する処理」、「森林領域の抽出」等を備える。
【０１６４】
（地上開度φ１´、地下開度φ２´の算出の見直し）
実施の形態１は疑似樹冠モデルという考えをもとに、DHM を用いて算出した地上・地下開度を樹冠モデルに合わせて数値を置き換えていた。このとき、樹冠端部から樹木頂点の角度を60°と設定していたが、本実施の形態２では現地で確認したスギ・ヒノキの樹形に合うように70°へと変更する。
【０１６５】
また、実施の形態１では林縁部の極大化を防ぐために、林縁部から地表を見下ろしたときに−60°以上の角度となる場合、林縁部の値を−60°に置き換えていた。しかし、依然として林縁部は凸の形状となり、第１−３の樹冠形状指数Ｇｅｉは大きく、一方、隣接する樹冠の境界に位置する場合は、凹の形状となり、第１−３の樹冠形状指数Ｇｅｉは小さくなる傾向にある。
【０１６６】
また、どちらも樹冠端部に位置するという点では一致しているものの、第１−３の樹冠形状指数Ｇｅｉに大きな違いが生じる。そこで、この違いを解消するため、新たな樹冠形状の−70°以上の角度となるときに林縁部と規定し、その値を70°(隣接する樹冠の境界に位置する場合の角度)に置き換える。なお、ここではDHM が3m 以下の非森林領域は樹冠形状指数の最小値（-70）とした。詳細については後述する。
【０１６７】
（局所加重平均）
実施の形態１では、算出された地上・地下開度（φ１、φ２）を用いて尾根谷度の算出式（φ3 = (φ1 −φ2 )/ 2）と同様な計算式で第１−１の樹冠形状指数φ3´ を求め、この第１−１の樹冠形状指数φ3´ を8 方位で計算して、平均化した第１−２の樹冠形状指数ｅｉを更に局所領域で平滑化処理して第１−３の樹冠形状指数Ｇｅｉとしていた。疎林において、この平滑化は複雑な樹冠形状を滑らかにして、同一樹冠内に生じる複数の樹木頂点を絞る効果があった。
【０１６８】
しかしながら、密林においては上記作用により樹木頂点の抽出数を減らすという結果になり、抽出精度が低くなる場合がある。隣接する樹木頂点が平滑化により樹木頂点数が減少する例を図２８に示す。図の上部分は樹冠形状指数の画像を数値化したものであり、下部分は横線部（２８ａ、２８ｂ）の断面形状を示している。
【０１６９】
図２８から、同一局所領域内に樹木頂点が複数含まれる場合、[１]平滑化により樹木頂点が減少する、[２]樹木頂点の位置座標が変わることが問題であるとわかる。このため、本実施の形態２では後述する局所加重平均によってこの問題を解決する。
【０１７０】
（森林領域の抽出）
樹冠形状指数は地表面（DSM）と地盤面（DEM）の差分（DHM）をもとにわずかな凹凸からも樹冠モデルを仮想的に作成するものである。一方で、樹木が存在しないギャップや林道などの領域においても樹冠モデルを作成してしまい、樹木が存在しない領域でも樹木頂点の抽出がなされてしまう。
【０１７１】
そこで、樹木が存在しない領域において樹木頂点の抽出を行わないように、森林領域の抽出を行った（図２９参照）。森林領域の抽出にはDHM を用い、下層植生を含まないようにここでは一例としてDHMが3m以上の領域を森林領域とした。
【０１７２】
図３０は本実施の形態２の樹木頂点認識装置の概略構成図である。図３１は本実施の形態２の樹木頂点認識装置を説明する概略フローチャートである。初めにこの概略フローチャートを説明する。
【０１７３】
図３１に示すように、航空レーザ計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ３０）。また、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ３１）。
【０１７４】
そして、DSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）を樹木高グリッドデータＤｄｉとして求めて（Ｓ３２）、樹冠形状の再現処理を行う（Ｓ３３）。
【０１７５】
この樹冠形状の再現処理は、樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）に対して樹木間隔程度（例えば１．５ｍ〜５ｍ範囲（図７を参照））を設定する（Ｓ３３ａ）。
【０１７６】
次に、各グリッドデータと自身のデータとのアスペクト比を求め、このアスペクト比を一定条件下で新たなアスペクト比（スギ、ヒノキの樹冠形状に合わせて７０°に変更）に合わせて変更し、変更したアスペクト比（（ΔＺｉ／ｄｉ）＝Ａｄｉ）を用いて第２の地上開度φ１´及び第２の地下開度φ２´を算出する（Ｓ３３ｂ）。
【０１７７】
本実施の形態２では現地で確認したスギ・ヒノキの樹形に合うように70°へと変更する。また、実施の形態１では、林縁部の極大化を防ぐために、林縁部から地表を見下ろしたときに−60°以上の角度となる場合、林縁部の値を−60°に置き換えていた。
【０１７８】
しかし、依然として林縁部は凸の形状となり、樹冠形状指数は大きく、一方、隣接する樹冠の境界に位置する場合は、凹の形状となり、樹冠形状指数は小さくなっていた。
【０１７９】
どちらも樹冠端部に位置するという点では一致しているものの、樹冠形状指数に大きな違いがあった。そこで、この違いを解消するため、新たな樹冠形状の−70°以上の角度となるときに林縁部と規定し、その値を70°(隣接する樹冠の境界に位置する場合の角度)に置き換える。
【０１８０】
そして、第２の地上開度φ１´及び第２の地下開度φ２´を尾根谷度の算出式により、第２−１の樹冠形状指数φ３´を算出してこれを８方位で平均化した第２−２の樹冠形状指数ｅｉを得た後で、ガウシアンフィルタを用いた平滑化（図３２参照）を行った第２−３の樹冠形状指数ｅｅｉを得る。次に、森林領域（例えばDHMが3m以上の領域）を規定し、この森林領域における第2-3の樹冠形状指数eeiはそのままの値とし、非森林領域（例えばDHMが3m未満の領域）の樹冠形状指数を任意の数値（例えば「-70」）とする第２−４の樹冠形状指数Ｅｅｉを求め（Ｓ３３ｃ）（図２９右図参照）、樹冠形状指数画像を得る（Ｓ３３ｄ）。
【０１８１】
この第２−４の樹冠形状指数Ｅｅｉ（単に樹冠形状指数Ｅｅｉという場合もある）と実施の形態１の第１−３の樹冠形状指数Ｇｅｉとを総称して樹冠形状指数ＥＧｉという。
【０１８２】
前述の開度の算出で、樹木間隔を樹木間隔程度としたが、これは検索範囲が長すぎると過小抽出となり、短すぎると過大抽出となるためである。一例として、表４に示すように立木密度から樹木間隔を求めるのが好ましい。
【表４】
【０１８３】
但し、Ｈ１〜Ｈ５は検証ポイント（所定の領域を有する）を示す。
【０１８４】
次に、樹木頂点の算出処理を行う（Ｓ３４）。
【０１８５】
この樹木頂点の算出処理は、第２−４の樹冠形状指数Ｅｅｉと動的閾値との比較により樹冠部のデータＤｈｉを抽出して（Ｓ３４ａ）、樹冠部画像を得る（Ｓ３４ｂ）。この樹冠部域の抽出を行う動的閾値の局所領域サイズは樹木間隔程度である。
【０１８６】
局所領域を樹木間隔程度とするのは、樹木毎の第２−４の樹冠形状指数Ｅｅｉに応じて樹冠の上部領域を抽出するためである。
【０１８７】
動的閾値は局所領域の平均値（μ）に標準偏差（σ）を足した値であり、この閾値に
より1 本の樹木の樹冠領域に対して約16％の中央部（樹冠の上部領域）を抽出することができる。
【０１８８】
また、樹冠部域の抽出は上記の通り局所領域で閾値を求め、この値と局所領域の中心画素の第２−４の樹冠形状指数Ｅｅｉとを比較して第２−４の樹冠形状指数Ｅｅｉが大きければ樹冠部域のデータＤｈｉ（ＤＱｉ）として抽出する。
【０１８９】
表9に樹木間隔と各プロットの局所領域サイズを示す。
【表９】
この局所領域サイズを用いて抽出した樹冠部域を図３４に示す。
【０１９０】
そして、樹冠部のデータのデータＤＱｉに対応する第２−４の樹冠形状指数Ｅｅｉの中から樹木冠形状指数が最高値となる第２−４の樹冠形状指数Ｅｅｉを検索し、これを樹木頂点Pipとする（Ｓ３４ｃ）。
次に、検索した樹木頂点の隣接程度を判断し、隣接していると判断したときは、第２−４樹冠形状指数Ｅｅｉの大小関係を基に隣接する樹木頂点のデータＰｉｐの第２−４の樹冠形状指数Ｅｅｉの大小関係をもとに、樹木頂点のデータＰｉｐを除去して（Ｓ３４ｄ）、樹木頂点Ｐｉとした画像（図３５参照）を得る（Ｓ３４ｅ）。
【０１９１】
次に、図３０の構成を説明する。
【０１９２】
（図３０の説明）
図３０においては図１と同一の符号のものについては説明を省略する。図３０に示すように、本実施の形態２の樹木頂点認識装置は、第２の樹冠形状指数算出部４１と、第２の樹木頂点決定部４３（隣接樹木除去処理を含む）等を備える。
【０１９３】
第２の樹冠形状指数算出部４１は、実施の形態１と同様な樹冠形状指数の算出処理（図示せず）を備えると共に、新たな第２のアスペクト比更新処理４１ａと、第２の樹冠形状指数計算処理４１ｂと、ガウシアンフィルタ処理４１ｃと、森林領域抽出処理４１ｄ等を備えている。
【０１９４】
第２のアスペクト比更新処理４１ａは、樹木高グリッドデータＤｄｉ（Ｘ、Ｙ、ＤＳＭとＤＥＭの差（ＤＳＭ））をファイル１５から全て読み込んで、これらのグリッドデータＤｄｉにおけるアスペクト比（樹木高Ｚｃｉ/距離ｄｉ）を求める。このアスペクト比は、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９５】
そして、図３７、図３８に示す様にアスペクト比を更新する。このアスペクト比が−70°以上の角度となるときに、その値を70°に置き換える。また、−７０°から０度の場合は−７０°に置き換える。この新たな更新アスペクト比は、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。また、本実施の形態２の新たなアスペクト比の更新については詳細に後述する。
【０１９６】
第２の樹冠形状指数計算処理４１ｂは、アスペクト比更新処理が終わると、３×３の領域毎にこの変更された新たな最大値、新たな最小値を用いて新たな第２の地上開度φ１´、第２の地下開度φ２´を求め、これらの開度から第２−１の樹冠形状指数φ３´を求めて、８方向で平均化して、第２−２の樹冠形状指数ｅｉを順次求める。
【０１９７】
ガウシアンフィルタ処理４１ｃは、第２−２の樹冠形状指数ｅｉの算出処理が終わると、３×３の領域毎の第２−２の樹冠形状指数ｅｉに対して局所加重平均化処理を行った第２−３の樹冠形状指数ｅｅｉを得る。この第２−３の樹冠形状指数ｅｅｉはグリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９８】
森林領域抽出処理４１ｄは、ガウシアンフィルタ処理後の第２−３の樹冠形状指数ｅｅｉの関連する樹木高グリッドデータのDHMが例えば３ｍ以上の場合に、森林領域の第２−４の樹冠形状指数Ｅｅｉとして抽出する。この第２−４の樹冠形状指数Ｅｅｉは、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９９】
閾値算出部３５は、表示部１３に表示された画像（グリッドデータ又はオルソフォト画像又はオルソフォト画像とグリッドデータとを重ね表示）内で、予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第２−４の樹冠形状指数Ｅｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【０２００】
樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の第２−４の樹冠形状指数Ｅｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・、図２２（ａ）参照）としてファイル２４に保存する。
【０２０１】
第２の樹木頂点決定部４３は、メモリ２４のデータＤｈｉの隣接するまとまり毎に番号を割り振ったデータＤＱｉ（樹冠形状指数のデータ）をメモリ２８に生成し、まとまり毎の番号の内で最大の第２−４の樹冠形状指数Ｅｅｉを持つ点を樹木頂点Pipとする。
【０２０２】
さらに、樹木頂点Pipに対して例えば3×3の局所領域を設定し、順次この局所領域を移動させ、局所領域の中央に樹木頂点Pipが位置するときに局所領域内に他の樹木頂点Pipが存在するかどうかを判定する。存在している場合は、第２−４の樹冠形状指数Ｅｅｉの大小を比較し、中央の樹木頂点Pipの樹冠形状指数Ｅｅｉが小さい場合、この樹木頂点を消去する。
【０２０３】
そして、隣接樹木消去処理が終わったデータを樹木頂点Ｐｉとしてメモリ２９に記憶する。
【０２０４】
出力部３０は、ファイル２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。
【０２０５】
上記のように構成された実施の形態の樹木頂点認識装置について以下にフローチャートを用いて説明する。
【０２０６】
本実施の形態２では、上記の実施の形態１と同様な樹冠立体角算出処理によって、既に検索方向ｋｉｐ毎のグリッドデータＤｐｐｉと樹木高グリッドデータＤｄｉｐとの更新アスペクト比Ａｄｉが求められてファイル３７に記憶され、アスペクト比更新処理が起動したとする。
【０２０７】
図３６は実施の形態２のアスペクト比更新処理を説明するフローチャートである。
【０２０８】
実施の形態２のアスペクト比の更新は図３７に示すように更新する。
【０２０９】
つまり、
（6）.「ΔZ/d」が-2.7474（-70 度）以下の場合は2.7474（70 度）とする。
【０２１０】
（樹冠端部の極大化防止）
（7）.「ΔZ/d」が-2.7474（-70 度）から0（0 度）の場合は-2.7474（-70 度）とする。
【０２１１】
（樹冠山地形の強調）
（8）.「ΔZ/d」が0（0 度）から2.7474（70 度）の場合は2.7474（70 度）とする。
【０２１２】
（樹冠谷地形の強調）
（9）.「ΔZ/d」が2.7474（70 度）以上の場合は「ΔZ/d」の値を据え置く。
【０２１３】
（樹木頂点の強調）
（10）.「ΔZ/d」が隣接する8 方位全てがマイナスである場合「ΔZ/d」を1000倍する。
【０２１４】
という具合に更新する。
【０２１５】
アスペクト比更新処理は、図３９に示すようにファイル３７の樹木高グリッドデータＤｄｉを設定し、検索方向ｋｉｐを設定する（Ｓ１００、Ｓ１０１）。
【０２１６】
次に、設定された番号の樹木高グリッドデータＤｄｉの検索方向ｋｉｐ内でのアスペクト比データViの中から値が最大なアスペクト比（ＭＡＸＡｄｉ）、値が最小のアスペクト比（ＭＩＮＡｄｉ）を抽出する（Ｓ１０２）。
【０２１７】
次に、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉが上記（８）の０〜２．７４７４の範囲（角度では０度〜７０度）かどうかを判定する（Ｓ１０３）。
【０２１８】
ステップＳ１０３において、０〜２．７４７４の範囲と判定したときは、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉを２．７４７４（ＭＡＸＡｄｉ、ＭＩＮＡｄｉ：Ａｄｂｉ）として処理をステップＳ１１１に移す（Ｓ１０４）。これにより、頂点から下の場合は、頂点から下の谷地形であることを強調する。
【０２１９】
また、ステップＳ１０３において、０〜２．７４７４の範囲（角度では０度〜７０度）ではないと判定したときは、読み込んだアスペクト比が上記（９）の２．７４７４以上かどうかを判定する（Ｓ１０５）。
【０２２０】
ステップＳ１０５において、アスペクト比が２．７４７４以上と判定したときはそのアスペクト比（Ａｄｉ）をそのままの値（Ａｄａｉ）にして処理をステップＳ１１１に移す（Ｓ１０６）。
【０２２１】
また、ステップＳ１０５において、読み込んだアスペクト比が２．７４７４以上ではないと判定したときは、読み込んだアスペクト比が上記の（６）の−２．７４７４以下かどうかを判定する（Ｓ１０７）。
【０２２２】
ステップＳ１０７において、−２．７４７４以下であると判定したときは、ステップＳ１０９に処理を移す。ステップＳ１０９においては、抽出したアスペクト比を＋２．７４７４（Ａｄbｉ）にする。
【０２２３】
これにより、樹冠端部と地盤との位置関係から生じるアスペクト比の極小化（樹冠形状指数では極大化）を防止する。
【０２２４】
また、ステップＳ１０７において、−２．７４７４以下ではないと判定したときは、抽出したアスペクト比（０〜−２．７４７４（７）(０を除く)）を−２．７４７４（Ａｄｃｉ）にする。これにより、樹冠端部から上の場合は、樹冠端部から上の山地形であることを強調する。
【０２２５】
また、前述のステップＳ１０４、ステップＳ１０６、ステップＳ１１０、ステップＳ１０９のアスペクト比をファイル３７に記憶する（Ｓ１１１）
次に、検索方向ｋｉｐが８かどうかを判定し（Ｓ１１２）、８方向になっていないときは検索方向を更新して処理をステップＳ１０１に戻す（Ｓ１１３）。
【０２２６】
ステップＳ１１３において、隣接する８方向の樹木高グリッドデータのアスペクト比が、水平軸を基準として全てマイナスと判定したときは、８方向の更新後のアスペクト比を１０００倍（Ａｄｄｉ）する（Ｓ１１５）。すなわち、頂点を強調する。
【０２２７】
次に、検索エリア内で樹木高グリッドデータＤｄｉが最後かどうかを判定し（Ｓ１１６）、最後でないときは次の番号の樹木高グリッドデータＤｄｉに更新して処理をステップＳ１００に戻す（Ｓ１１７）。
【０２２８】
すなわち、ファイル３７には、ステップＳ１０９、Ｓ１１０の処理によって、図３８に示すように、設定された番号の樹木高グリッドデータＤｄｉと検索方向ｋｐｉと抽出したＭＡＸＡｄｉ、ＭＩＮＡｄｉと更新アスペクト比とが対応させられて保存されている。
【０２２９】
更新アスペクト比は、アスペクト比Ａｄｉのままの値（Ａｄａｉ）又は２．７４７４（Ａｄｂｉ）又は−２．７４７４（Ａｄｃｉ）若しくは−２．７４７４を１０００倍したアスペクト比（Ａｄｄｉ）のいずれかである。
【０２３０】
従って、第２の樹冠形状指数計算処理４１ｂは、Ｄｄｉ、ｋｉｐに対応する更新アスペクト比を引当、この更新アスペクト比の中から最大アスペクト比ＭＡＸＡｄｉ（アスペクト比Ａｄｉのまま又は２．７４７４又は−２．７４７４又は−２７４７．４のいずれか）と、最小アスペクト比ＭＩＮＡｄｉ（アスペクト比Ａｄｉのまま又は２．７４７４又は−２．７４７４又は−２７４７．４のいずれか）読み込むことになる。
【０２３１】
（第２の樹冠形状指数計算処理）
そして、第２の樹冠形状指数計算処理４１ｂは、第２の地上開度φ１´と第２の地下開度φ２´を算出する。
【０２３２】
第２の地上開度φ１´＝９０度−tan-1(ＭＡＸＡｄｉ)
第２の地下開度φ２´＝９０度＋tan-1(ＭＩＮＡｄｉ)
次に、この第２の地上開度φ１´と第２の地下開度φ２´とを用いて第２−１の樹冠形状指数φ３´を求める。
【０２３３】
第２−１の樹冠形状指数φ３´＝（φ１´−φ２´）／２
次に、８方位全てに対して、これを平均化（８方位）し、この平均化された第２−１の樹冠形状指数φ３´を、樹木高グリッドデータＤｄｉを中心とした検索範囲での第２−２の樹冠形状指数ｅｉとして、樹冠グリッドデータＤｄｉにリンク付けしてファイル３７に記憶する。
【０２３４】
最後に第２−２の樹冠形状指ｅｉを検索範囲に対応させてファイル３７に記憶する。
【０２３５】
次に、本実施の形態２で用いるガウシアンフィルタ処理４１ｃについて説明する。
【０２３６】
この手法は局所平均（これまでの平滑化）と異なり、対象局所領域に重みを持たせ、画像をより鮮明にすることができる。局所平均と局所加重平均の例を図４０に示す。図４０（ａ）が局所平均であり、図４０（ｂ）が局所加重平均の処理結果である。
【０２３７】
つまり、本実施の形態ではガウシアンフィルタという局所加重平均を用いて第２−２の樹冠形状指数ｅｉを平滑化する。
【０２３８】
ガウシアンフィルタは下記の式（E）のガウス関数で加重を決定する。
【数１】
【０２３９】
ここで、G は各画素の加重、i,j は縦横それぞれの中心画素からの距離（画素）、νは変数を表す。
【０２４０】
ガウシアンフィルタは変数νを任意で変えることで、画像の平滑化の度合いを変化させることができるという利点を持っている。
【０２４１】
一般的に局所領域の大きさが(2N＋1)×(2N＋1)
の場合、νはN/2 とするのが目安である。また、νの値が小さいとより原画像に近く、νの値が大きいと局所平均画像に近づく特徴がある。図４４にガウシアンフィルタを用いた平滑化の例を示す。なお、νは目安となっているN/2 とした。図から上記で述べた局所平均による平滑化の2つの問題点を解決できていることから、本実施の形態２ではガウシアンフィルタを用いた平滑化を採用する。νは、０．５〜０．７（図３２参照）とする（好ましくは０．７）。
【０２４２】
ガウシアンフィルタ処理４１は、このような第２−３の樹冠形状指数ｅｅｉを、メモリ３７に記憶し、森林領域抽出処理４１ｄが以下の処理を行う。
【０２４３】
（森林領域抽出処理）
次に、森林領域抽出処理を説明する。
【０２４４】
樹冠形状指数は地表面（DSM）と地盤面（DEM）の差分（DHM）をもとにわずかな凹凸からも樹冠モデルを仮想的に作成するものである。一方で、樹木が存在しないギャップや林道などの領域においても樹冠モデルを作成してしまい、樹木が存在しない領域でも樹木頂点の抽出がなされてしまう。
【０２４５】
そこで、樹木が存在しない領域において樹木頂点の抽出を行わないように、森林領域の抽出を行った（図２９参照）。森林領域の抽出にはDHM を用い、下層植生を含まないようにここでは一例としてDHMが3m以上の領域を森林領域とした。そして、森林領域の第2-3の樹冠形状指数eeiはそのままの値を用い、非森林領域（例えばDHMが3m未満の領域）の樹冠形状指数を任意の数値（例えば「-70」）にしたものを第２−４の樹冠形状指数Ｅｅｉとした。この第２−４の樹冠形状指数Ｅｅｉは、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０２４６】
このため、ファイル３７には図３９に示すように、設定された樹木高グリッドデータＤｄｉと、検索方向ｋｉｐと、更新前のＭＡＸ、ＭＩＮのアスペクト比と、更新後のＭＡＸ、ＭＩＮのアスペクト比とが検索方向毎に、リンクつけされて記憶されている。
【０２４７】
また、ファイル３７には図３９に示すように、設定された樹木高グリッドデータＤｄｉの番号と、そのＸ、Ｙ座標と、樹木高の差Ｚと、次の番号の樹木高グリッドデータと、このＸ、Ｙ座標と、距離ｄｉと、第２―１の樹冠形状指数φ３´と、この第２―１の樹冠形状指数φ３´を得たときの第２の地上開度φ１´と、第２の地下開度φ２´とが検索方向にリンクつけされ、また、これらの８方向全ての平均化された第２−２の樹冠形状指数ｅｉ及びガウシアンフィルタ後の第２−３の樹冠形状指数ｅｅｉ、森林領域抽出後の第２−４の樹冠形状指数Ｅｅｉ等が対応させて記憶されている。これらを総称して第２の樹冠形状指数情報Ｇｉｉともいう。
【０２４８】
次に、閾値算出部３５が予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第２−４の樹冠形状指数Ｅｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【０２４９】
次に、樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の第２−４の樹冠形状指数Ｅｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル２４に保存する。
【０２５０】
次に、第２の樹木頂点決定部４３は、メモリ２４のデータＤｈｉの隣接するまとまり毎に番号を割り振ったデータＤＱｉ（樹幹形状指数のデータ）をメモリ２８に生成する（図４６（ａ））。
【０２５１】
このデータＤＱｉに対して以下の樹木頂点抽出処理を行う。
【０２５２】
データＤＱｉは樹冠領域毎（まとまり毎）に重複しない番号が割り振られる。一つの番号につき、複数のデータDdiが存在する場合、この番号の内で最大の樹冠形状指数を持つ点を樹木頂点Pipとする。
【０２５３】
さらに、樹木頂点が隣接し、過剰抽出となる場合があるため、隣接樹木頂点の除去処理を行う。
【０２５４】
樹木頂点抽出結果（上記手法により抽出された樹木頂点Pipの抽出精度）からスギ林だけでなく、ヒノキ林においても精度良く抽出する手法が開発できたと考える。なお、樹冠部域画像と樹木頂点抽出結果を参照すると、抽出樹木頂点が隣接している場合に若干多めに樹木頂点を抽出する傾向にあるが、この点に対しては以下の手法により解決できる。
【０２５５】
図３５に1 つの樹冠に複数の樹木頂点が存在している状況を示す。これは、樹冠部域が1 つの樹冠に対して1 つの樹冠部域が抽出されず、樹冠部が分割されたために生じたと考えられる0.5m グリッドにおいて樹木頂点が斜めに隣接する場合、2 点間の距離は0.7m であり、立木密度を考慮するとその間隔は短い。このことからも1 つの樹冠から過剰に樹木頂点が抽出されていることがわかる。
【０２５６】
この問題を解決するために、樹木頂点の抽出でも使用した樹冠形状指数の大小関係を利用して、隣接する樹木頂点を除去した。具体的には、Pipに対して例えば3×3の局所領域を設定し、順次この局所領域を移動させ、局所領域の中央に樹木頂点Pipが位置するときに局所領域内に他の樹木頂点Pipが存在するかどうかを判定し、存在している場合は、樹冠形状指数Ｅｅｉの大小を比較し、中央の樹木頂点Pipの樹冠形状指数Ｅｅｉが小さい場合、この樹木頂点を消去する。
【０２５７】
図４５は隣接樹木除去処理を説明するフローチャートである。
【０２５８】
隣接樹木処理は、樹木頂点のデータPipに対して３×３のフィルタをかけて、このフィルタの中央の位置ＪＤｉを設定する（Ｓ２０１、図４６（ｂ））。
【０２５９】
次に、中央の位置ＪＤｉの隣接8方向の領域ＪＤｒｉに樹木頂点のデータPipがあるかどうかを判定する（Ｓ２０２）。
【０２６０】
ステップＳ２０２において、データＰｉｐがあると判定したときは、隣接8方向の領域ＪＤｒｉにデータＰｉｐに対応する第２−４の樹冠形状指数Ｅｅｉと中央の位置ＪＤｉの第２−４の樹冠形状指数Ｅｅｉとを比較して値の大小を判定する（Ｓ２０３、Ｓ２０４）。
【０２６１】
ステップＳ２０４において、ＪＤｉの樹冠形状指数Ｅｅｉが小さいと判定したときは、ＪＤｉの樹木頂点のデータPipを「0」とする（Ｓ２０６）。
【０２６２】
また、ステップＳ２０４において、ＪＤｉの樹冠形状指数Ｅｅｉが大きいと判定したときは、樹木頂点のデータPipがメモリ２８に他にあるかどうかを判定する（Ｓ２０７）。他にあると判定したときは、ＪＤｉを次の番号に更新して処理をステップＳ２０１に移す（Ｓ２０９）。
【０２６３】
図４６（ｂ）においては、３×３フィルタの中央に「２」が存在し、斜め方向に「１」が存在している。「２」に対応する第２−４の樹冠形状指数が、「１」に対応する第２−４の樹冠形状指数より小さい場合は、図４６（ｃ）に示すように「２」は消去される。
【０２６４】
また、ステップＳ２０７において、データＰｉｐが他にないと判定したときは、メモリ２８の全ての樹木頂点のデータPipを樹木頂点Ｐｉとしてメモリ２９に記憶する（Ｓ２０８）。
【０２６５】
次に、樹木頂点消去について説明を補充する。
図４４ に樹木頂点を除去した結果を示す。図に示す通り、同一樹冠内の樹木頂点が除去され、1つの樹冠に対して、1 つの樹木頂点が抽出されることが確認できた。
【０２６６】
隣接する樹木頂点を除去した結果を表５に示す。
【表５】
【０２６７】
H4-1、H4-2、H5-1 以外のプロットでは抽出精度に違いはないが、特にH4-1、H4-2 では樹頂点の抽出精度が向上した。現地で確認した樹冠形状とオルソ画像を勘案すると、H4 プロットは樹冠表層が他のプロットに比べて凹凸が明瞭ではなく、複雑な形状をしていると思われる。このため、1 つの樹冠に対して隣接する２つの樹冠上部が抽出され、その結果、１つの樹冠に対して隣接する２つの樹頂点が抽出されてしまったと考えられる。
【０２６８】
＜実施の形態３＞
前述の実施の形態１及び実施の形態２のＤＳＭは、 樹冠形状明確化を行うために樹冠内部の標高データを除去している。
【０２６９】
これは、樹冠を透過するパルスや航空機から斜めにレーザーが照射されることで樹木間の隙間から地盤に到達するパルスが存在するため、従来のDSM には樹冠内部の標高データ（樹冠表層の標高データではなく、樹冠等を透過して取得された標高データ）も含まれている（図４７参照）。このDSM を用いて樹高データ(DHM)を算出すると、樹冠形状が不明瞭になってしまう傾向にある。
【０２７０】
従って、樹木頂点抽出精度を向上させるためには、樹冠内部の標高データを除去し、樹冠形状を明確化したDSM を作成する必要がある。
【０２７１】
手動により樹冠内部の標高データを除去し、樹冠形状を明確化する手法もある。しかし、広範囲にわたる森林でのデータの処理には時間と労力がかかり、かつ、人為的な要因による再現性の問題が生じる可能性がある。
【０２７２】
先述したとおり、従来のDSM では地盤面の標高データも含まれているため、計測データは図４８（ａ）のようになっている。赤点で示す樹冠表層の標高データと青点で示す地盤面の標高データではその分布域に違いがあるので、この差を用いて地盤面に存在する標高データを除去することができる。
【０２７３】
反射強度は反射する物体により、その値が異なり、樹木などの植物の反射強度は大きいが、道路などの人工物や地盤などの土の反射強度は小さいという特徴を持つ。また、樹冠などを透過することで、反射強度は減衰する。この特徴を活かし、標高データだけでは除去できない樹冠内部に存在する標高データを除去する。
【０２７４】
図４８（ｂ）に示す青点は樹冠内部に存在する標高データであり、標高データだけでは除去できないデータである。赤点で示す標高データは樹冠表層に存在するため、反射強度は大きい。一方、青点で示す樹冠内部の標高データは樹冠を透過し、減衰しているために反射強度は小さくなる。この反射強度の違いを用いて樹冠内部のデータを除去することができる。
【０２７５】
以上より、自動処理による樹冠形状明確化に十分な成果が得られたので、本手法による樹冠形状明確化DSM を作成し、以後の解析（樹木頂点の抽出等）を行うのが好ましい。
【０２７６】
なお、本実施の形態１および実施の形態２は下記に説明するように利用されることもある。
【０２７７】
（一元管理化）
本実施の形態２により高精度に林分パラメータが取得できることがわかった。
【０２７８】
任意の林分に対して、「平均樹高」「立木密度」「推定樹冠面積」といった林分パラメータが得られる。さらにこのデータを用いることで間伐林の選定に役立つ「収量比数」や「適正立木密度」が得られる。また、レーザー計測データからは「林分傾斜角」や「レーザーパルスの地盤到達率」といった情報も得られる。データの一元管理化した例を以下に示す。
【０２７９】
[１］林分に対して複数の情報を付与できる。
【０２８０】
[２］上記の取得情報を取得情報の位置や地図等の情報と共に可視的に管理できる。
【０２８１】
[３]データ管理・更新が容易で、データの信頼性が高められる。
【０２８２】
[４]一元データ管理により、データの比較・解析・評価が容易である。
【０２８３】
[５]他の情報（森林簿）と併せてデータを集約し、その情報も一元管理できる。
【０２８４】
以下に収量比数、森林整備ランクについて述べる。
【０２８５】
1）収量比数
収量比数は密度管理図をもとに算出するもので、最多密度における幹材積に対する現存幹材積の割合である。具体的には樹高と立木密度から算出される。梨のスギ、ヒノキ林の密度管理図に描かれている曲線式を次に示す。
【０２８６】
スギ林
V = (0.07155968H-1.373859 + 5062H-2.869785 / N )-1 式（F）
logNRF = 5.370947 −1.495926logH 式（G）
ヒノキ林
V = (0.035147H-1.080773 + 4711.2H-2.922894 / N )-1 式（H）
logNRF = 5.7384 −1.8482121logH 式（I）
NRF：最多密度におけるha あたりの本数（本/ha）である。
【０２８７】
本研究ではプロット内に含まれる樹木本数（現地調査結果、樹冠未到達の中下層木も含む）と抽出本数（開発手法によるもので、表層木のみとなる）とをもとにした立木密度を用いて収量比数を算出した。その結果を図４９に示す。この図からプロット間の相対的な大小関係は大きく変わらないため、間伐の優先度を選定する情報として表層木のみの収量比数を用いることが可能とわかる。
【０２８８】
また、収量比数を0.7 として適正立木密度を算出した。この数値は「中庸仕立
て」と呼ばれる一般的な間伐の指標である。
【０２８９】
森林整備ランク
森林の生育状況だけでなく、樹種や森林の生育している状況にも注目して森林整備の優先度を評価する指標が森林整備ランクである。これは「森林管理ランク（収量比数）」「森林特性ランク（樹種や林齢）」「国土保全ランク（傾斜や保全対象）」を求め、総合的に森林整備の優先度を把握することができる。
【０２９０】
森林整備ランクを算出した結果を図５０，５１に示す。このことから、スギ林においては本指標により劣勢木の潜在状況も加味した森林整備計画が可能と考える。
【０２９１】
また、取得した情報を林分毎に整備し、GIS データとして一元管理する例を図５２ に示す。
【０２９２】
これにより、森林整備ランクをオルソ画像情報に載せて把握したり、プロットを選択することで図中右上のような林分の詳細データを把握したりできる。成可能であり、下記の様な利用が考えられる。
【０２９３】
[１] 森林の状態や森林所有者の把握が容易
[２] 優先度を考慮した効率的な森林整備計画
[３] 情報共有・公開社会に適した森林情報のデータ整備・運用
詳細な現地調査を実施することなく、簡易・広範囲
に森林整備の優先度が把握できることは、林業経営の観点から求められているコスト削減や要間伐林の選定が急務となっている現在の社会情勢に貢献するものと考える。また、林分パラメータだけでなく、レーザー計測データのオルソ画像や既存森林簿等も情報として加えることで、さらにGIS 一元管理の利便性が高まり、公開・共有情報化社会に貢献するものと考える。
【０２９４】
さらに、災害等により、地面が削りとられた場合に、本手法により予めこの個所の付近のそう本数を得ておいて、再びこの個所の付近の本数を求めて、その差を求めると、どの程度の木が川に流れたかが分る。
【０２９５】
（樹冠面積の推定）
樹冠面積推定は樹冠形状指数Ｅｅｉをもとに、watershed アルゴリズムで行っている。実施の形態１の樹冠面積推定手法では、あらかじめ樹冠サイズを設定し、樹冠面積を推定していたため、疎密度にばらつきのある林分においては十分な樹冠面積の推定ができていなかった。そこで、アルゴリズムの見直しを行い、推定精度の向上を試みた。
【０２９６】
樹冠面積推定手法の改良は実施の形態１と同様にwatershed アルゴリズムを用いて行った。実施の形態１では樹木頂点から周囲を検索して集水域を推定していたが、本利用例では閾値以上の任意の点がどの樹木頂点の集水域に属しているかという考え方に基づいて樹冠面積を推定している。
【０２９７】
本利用例では、樹冠形状指数Eeiに設定する閾値以上を対象に樹冠面積を推定するため、最適な閾値を設定する必要がある。そこで、最適な閾値を設定するために、樹高データと樹冠形状指数の断面形状をもとに最適値を決める。
【０２９８】
樹高（DHM）3m 以上を森林領域と定義したことから、これに基づき、樹高3m 以上の領域と樹冠形状指数のある値以上の領域とがほぼ合致するとき、この樹冠形状指数Ｅｅｉの値を閾値とする。図３３に樹高と樹冠形状指数の2 値化画像（背景にオルソ画像を表示）を示している。
【０２９９】
図３３で示したように樹高3m以上の領域に相当する樹冠形状指数は-60 であったため、樹冠形状指数Eeiの閾値を-60 として樹冠面積を推定する。
【０３００】
樹冠面積推定結果と精度改良した手法を用いて推定した樹冠面積を実施の形態１で作成した推定樹冠面積と比較した結果を図４２に示す。図を見ると実施の形態１の推定面積に比べて、樹冠の大きさに応じた推定面積となっている。また、目視判読の樹冠面積と同様の大きさであり、推定手法が改善されていることがわかる。
【０３０１】
本利用例をヒノキ林に適用した結果の画像を図４３に示し、検証結果を表６、表７、表８に示す。
【表６】
【表７】
【表８】
【０３０２】
表６から樹木の抽出数の精度は目視判読と比較してH4 林分で十分な精度ではなく（密林では目視判読結果は過小抽出となる）、また、被覆面積率（樹冠面積合計の抽出割合）は目視判読に比べて、全体的に過大であった。しかし、調査林分の状況やfirst パルスの地盤到達率からヒノキ林分は鬱閉していると考えられ、目視判読の樹冠面積が過小と考えられる（目視判読の樹冠面積は樹冠の大きさ、ヒノキの樹冠形状、判読写真の解像度の問題から、鬱閉状態を表すような判読面積に上方修正することは困難）。このため、林分面積をもとに樹冠面積の抽出精度を検証した。その結果、林分面積に対して推定樹冠面積の抽出率は概ね100%となり、推定樹冠が鬱閉状態を表していることがわかる。
【０３０３】
表７には単樹冠面積の比較を示したが、密林では目視判読が過小抽出となることや推定樹冠面積が目視判読に比べ過大であることから十分な精度が得られていない。
【０３０４】
また、目視判読の抽出精度を利用して林分内に含まれる実樹木本数を推定した結果及び、推定実樹木本数と開発手法による抽出樹木本数との比較結果を表８に示す。目視判読の抽出精度から推定する実樹木本数と開発手法による抽出樹木本数との割合は約9〜11 割となり、林分全域で高精度に樹木本数が抽出できていると考えられる。
【０３０５】
さらに、このことと推定樹冠が鬱閉状態を再現できていたことを考えると表７ の推定単樹冠面積の再現性は高い。
【図面の簡単な説明】
【０３０６】
【図１】本発明にいたる経緯を説明の説明図である。
【図２】本発明にいたる経緯を説明する説明図である。
【図３】本発明にいたる経緯を説明する説明図である。
【図４】本実施の形態による尾根谷画像と樹冠形状指数画像との比較説明図である。
【図５】本実施の形態の樹木頂点認識装置の概略構成図である。
【図６】本実施の形態を説明する概略フローチャートである。
【図７】検索範囲の設定を説明する説明図である。
【図８】本実施の形態の樹冠形状指数による画像を説明する説明図である。
【図９】樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【図１０】樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【図１１】検索範囲における樹木高グリッドデータの抽出を説明する説明図である。
【図１２】更新するアスペクト比と角度の関係を説明する説明図である。
【図１３】アスペクト比の更新処理を説明するフローチャートである。
【図１４】抽出したＤｐｐｉと更新したアスペクト比等の情報を説明する説明図である。
【図１５】ファイル３７の樹冠形状指数情報Ｇｉの説明図である。
【図１６】樹冠形状指数の算出のフローチャートである。
【図１７】樹冠領域に対しての樹冠部との関係を説明する説明図である。
【図１８】局所領域の平均値（μ）と標準偏差（σ）との関係を説明する説明図である。
【図１９】具体的な樹冠部の抽出処理を説明するフローチャートである。
【図２０】頂点候補Ｄｇｉの説明図である。
【図２１】樹冠の上部のデータＤｈｉの説明図である。
【図２２】樹冠部の説明図である。
【図２３】本実施の形態による疎林の頂点抽出結果を説明する説明図である。
【図２４】本実施の形態の樹冠形状指数画像と尾根谷度画像との違いを説明するための説明図である。
【図２５】本実施の形態の樹冠形状指数画像と尾根谷度画像との違いを説明するための説明図である。
【図２６】本実施の形態による疎林の樹冠部抽出結果を説明する説明図である。
【図２７】本実施の形態の具体的なハードウエア構成図である。
【図２８】実施の形態２における樹冠形状指数の平滑処理（局所領域）による樹木頂点の減少を説明する説明図である。
【図２９】実施の形態２による森林領域抽出結果を説明する説明図である。
【図３０】本実施の形態２の樹木頂点認識装置の概略構成図である。
【図３１】本実施の形態２を説明する概略フローチャートである。
【図３２】本実施の形態２のガウシアンフィルタと局所平均の樹木頂点抽出結果の比較を説明する説明図である。
【図３３】本実施の形態２の利用例（樹冠面積の推定）の森林領域の最適閾値の設定根拠を説明する説明図である。
【図３４】本実施の形態２による樹冠部域の抽出結果を説明する説明図である。
【図３５】隣接頂点の不具合を説明する説明図である。
【図３６】本実施の形態２のアスペクト比更新処理を説明するフローチャートである。
【図３７】本実施の形態２の地上・地下開度の置き換え変更を説明する説明図である。
【図３８】本実施の形態２の抽出したＤｐｐｉと更新したアスペクト比等の情報を説明する説明図である。
【図３９】実施の形態２のファイル３７の樹冠形状指数情報Ｇｉｉの説明図である
【図４０】局所平均と局所加重平均を説明する説明図である。
【図４１】本実施の形態２のガウシアンフィルタの説明図である。
【図４２】本実施の形態２の利用例（樹冠面積の推定）の樹冠面積推定結果を説明する説明図である。
【図４３】本実施の形態２の利用例（樹冠面積の推定）の樹冠面積推定結果を説明する説明図である。
【図４４】本実施の形態２の隣接頂点樹木前後の比較を説明する説明図である。
【図４５】本実施の形態２の隣接樹木頂点除去処理を説明するフローチャートである。
【図４６】本実施の形態２の隣接樹木頂点の除去を説明する説明図である。
【図４７】実施の形態３のＤＳＭの標高データの削除を説明する説明図である。
【図４８】実施の形態３の樹冠内部データを説明する説明図である。
【図４９】本実施の形態２の利用例（収量比較）を説明する説明図である。
【図５０】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【図５１】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【図５２】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【符号の説明】
【０３０７】
１０ データベース
１１ データベース
１３ 表示部
１４ 樹木高算出部
１７ 尾根谷度算出部
２１ 樹木頂点候補抽出部
２３ 樹冠部抽出部
２５ 第１の樹木頂点決定部
３０ 出力部
３２ 第１の樹冠形状指数算出部
３５ 閾値算出部
４１ 第２の樹冠形状指数算出部
４１ａ 第２のアスペクト比更新処理
４１ｂ 第２の樹冠形状計算処理
４１ｃ ガウシアンフィルタ
４１ｄ 森林領域抽出処理
４３ 第２の樹木頂点決定部
【技術分野】
【０００１】
本発明は、レーザ計測で得た多数のデータ（ＤＳＭ、ＤＥＭ）から単木の頂点を認識する樹木頂点認識装置に関する。
【背景技術】
【０００２】
森林のCO2吸収、水源涵養、木材生産などの様々な機能を向上させるためには森林管理が必要である。森林管理を行うに当たっては、樹木本数などの林分パラメータが必要となる。
【０００３】
例えば、山地の特定エリアの唐松の本数を特定する場合には、写真画像等のそのエリア内の唐松の樹木の頂点を目視で特定し、これらの総数を求めていた。
【０００４】
また、レーザ測距儀（以下、航空レーザ）を活用して樹木本数を推定する場合もある。
【０００５】
具体的には、航空レーザ計測により取得されたデータを用いて表層のグリッドデータ（以下、ＤＳＭ：Digital Surface Model）と地盤のグリッドデータ（以下、ＤＥＭ：Digital Elevation Model）を作成する。さらに、DSMとDEMデータの差分をとることによって樹木高のグリッドデータを算出し、樹木高や樹木高の極大値を判別指標として単木抽出を行っていた。
【０００６】
また、特開２００４−３３１４９号公報（特許文献１）は、施業林分における立木の直径、本数及び配置の調査を自動的に行うものである。
【特許文献１】特開２００４−３３１４９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、写真画像から目視によって樹木本数を求めた場合は、非常に時間がかかる。特に、エリアが広くなるとそれだけ時間がかかる。また、写真画像から樹木を１本毎に目視で正しく特定していくには、経験が必要である。
【０００８】
また、従来のＤＥＭ、ＤＳＭによって単木抽出を行う手法は、樹冠形状が複雑であることから、一つの樹冠内には複数の極大値が存在しやすく、従来の手法で単木抽出処理を行うと実本数に比べ過剰に抽出されるデメリットがある。
【０００９】
本発明は以上の課題を解決するためになされたもので、自動的に単木の頂点を正しく認識する樹木認識装置を得ることを目的とする。
【課題を解決するための手段】
【００１０】
本発明の樹木頂点認識方法は、上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識方法である。
【００１１】
コンピュータが、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶手段に記憶するステップと、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出ステップと、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出ステップと、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出ステップと、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識ステップと
を行うことを要旨とする。
【発明の効果】
【００１２】
以上のように本発明によれば、樹冠形状指数によって、自動的に精度よく単木の頂点を認識できる。
【００１３】
このため、広い森林等を人手によって樹木一本毎に検索しなくともよいので、例えば森林の財産管理等に適用できる。
【発明を実施するための最良の形態】
【００１４】
＜実施の形態１＞
初めに経緯を説明する。
【００１５】
尾根谷度の算出は、次の様に行う。樹木高のグリッドデータを用いて、設定した対象範囲内（検索範囲）における天頂から地表（樹冠含む）に接するまでの角度と鉛直下方から地表に接するまでの角度を求める。そして、これらの角度を例えば、8方向で計算し、それぞれ平均する。その結果をそれぞれ地上開度（φ1）と地下開度（φ2）とする。地上開度と地下開度にはさまれる角度の2 等分線と水平線との角度を尾根谷度（φ3）とする（図１参照）。検索範囲は1m・2m・5m・10mの4 通り試算し、樹冠形状を適切に表現できた2m に決定した。
【００１６】
尾根谷度は大まかな樹冠形状を表現するものであり、樹高に比べ、樹木頂点と樹冠端部の値の差が大きく、樹木頂点において値が同じような傾向を示す（図２参照）。これにより画一的な閾値を用いて、樹冠部域の抽出が可能で、ここに尾根谷度の有用性があると考える。
【００１７】
次に、尾根谷度に局所最大値フィルタを適用し、局所領域内の最大値を中央画素値に置き換える。この際、局所領域を3×3 に設定した。このフィルタ処理後の尾根谷度と処理前の尾根谷度を比較し、処理前後で尾根谷度が同値であったときに、その地点を「樹頂候補点」として抽出する。
【００１８】
一方で、「尾根谷度＞0」は尾根地形を表現し、「尾根谷度＜0」は谷地形を表現するといった尾根谷度の性質を利用して、尾根谷度データに森林の状態に応じた閾値を設けることによって樹冠部を抽出する。そして、樹頂候補点と樹冠部領域を重ね合わせて、各樹冠部領域内において樹高が最も高い樹頂候補点を樹木頂点として抽出する。
【００１９】
しかしながら、抽出された樹木頂点の総数を現地の立木本数として検証していたが、抽出された樹木頂点は実際の立木本数よりも過小になる場合があることがわかった。すなわち、これまでの手法では抽出本数が過小になる可能性がある。
【００２０】
例えば、図３に示すように、尾根谷度は樹高データをもとに計算した地上開度と地下開度から算出される。このため、樹冠が隣木に接していない樹木は樹冠の端でDHM が急激に変化し、尾根谷度が大きくなる。これにより、樹冠端部が樹木頂点より大きな尾根谷度になる。
【００２１】
図４は尾根谷度画像と本実施の形態２の樹冠形状指数との比較を説明する説明図であり、図４（ａ）は尾根谷度が樹冠端部において極大化（画像が明るい領域）している例を示す画像である。
【００２２】
この極大化により、林縁部では、樹冠端部が頂点として認識されてしまったり、樹木頂点が樹冠内に複数生じたりする可能性がある。この問題を解決するには、樹冠端部の極大化を防ぎ、樹木の頂点の尾根谷度を強調して隣木との区別を明確にする方法が好ましい。
【００２３】
図５は実施の形態の樹木頂点認識装置の概略構成図である。図６は本実施の形態を説明する概略フローチャートである。実施の形態の樹木頂点認識装置の概略構成図を説明する前に、図６の概略フローを説明する。図６に示すように、航空レーザ計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ２０）。また、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ２１）。
【００２４】
そして、DSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）を樹木高グリッドデータＤｄｉとして求めて（Ｓ２２）、樹冠形状の再現を行う（Ｓ２３）。
【００２５】
樹冠形状の再現は、この樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）を用いて例えば５ｍ範囲（図７を参照）における各グリッドデータと自身のデータとのアスペクト比を求め、このアスペクト比を一定条件下で新たなアスペクト比に変更し、変更したアスペクト比を用いて地上開度φ１´及び地下開φ２´度を算出し（Ｓ２３ａ）、図４（ｂ）に示す樹冠形状指数画像のデータを算出する（Ｓ２３ｂ）。この樹冠形状指数画像については後述する。
【００２６】
そして、この樹冠形状指数Gｅiから樹木の頂点を抽出する（Ｓ２４）。樹木頂点の抽出は、樹冠形状指数Gｅiに局所最大値フィルタ（3×3範囲）を適用して（Ｓ２４ａ）、「樹木頂点の候補点」を抽出する（Ｓ２４ｂ）。
【００２７】
また、本実施の形態では、後述する動的閾値により樹冠部を抽出して（Ｓ２４ｃ）、樹冠部画像を得る（Ｓ２４ｄ）。
【００２８】
そして、頂点候補が存在する樹冠部領域毎に樹木頂点候補の中から最高の樹木高を有する点を求め、これを樹木頂点とし（Ｓ２４ｅ）、抽出する（Ｓ２４ｆ）。これによって得られた画像を図８に示す。
【００２９】
次に、図５の構成を説明する。
【００３０】
（図５の説明）
図５においては、表示処理部１９は、表示部１３（画面）に各データベース、メモリ等の画像を表示させたり、画面の画像データを読み込む等の機能を有しているが、各データベース、メモリ等とのやり取りを示す線については記載を省力している。
【００３１】
また、図５においては、特定地域の地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を記憶したデータベース１０と、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model）を記憶したデータベース１１と、特定地域のデジタル写真（オルソフォト画像）を記憶したデータベース１２と、樹木高算出部１４と、樹木高データが記憶されるファイル（メモリ）１５と等を備えている。
【００３２】
また、検索範囲決定部１７と、表示処理部１９と、局所最大フィルタ部２０と、樹木頂点候補抽出部２１と、樹木頂点候補データが記憶されるファイル２２と、樹冠部抽出部２３と、樹冠部のデータが記憶されるファイル２４と、第１の樹木頂点決定部２５と、樹木頂点データが記憶されるファイル２９と出力部３０等を備えている。
【００３３】
さらに、本実施の形態の樹木頂点認識装置は、第１の樹冠形状指数算出部３２と、閾値算出部３５と、局所領域設定部３６と等を備えている。
【００３４】
データベース１０は、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model：本実施の形態では０．５ｍ間隔の格子状のデータ）を記憶している。このグリッドデータは、地盤のＸ、Ｙ、Ｚ（地盤高Ｚａ）とを対応させて記憶させられている。つまり、レーザデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて地面を復元したものである。
【００３５】
データベース１１は、特定地域の表層（森林）のグリッドデータ（ＤＳＭ：Digital Surface Model：本実施の形態では０．５ｍ間隔）を記憶している。このグリッドデータ（本実施の形態では０．５ｍ間隔の格子状のデータ）は、樹木表面のＸ、Ｙ、Ｚ（Ｚｂ：樹木表面高）とを対応させて記憶している。つまり、レーザデータを読み込み、それぞれの同じ標高値を結んだ等高線図に対してＴＩＮを発生させて表層を復元したものである。
【００３６】
すなわち、これらのグリッドデータは、格子間隔０．５ｍのＸ−Ｙ平面直角座標系に、緯度経度に対応するＸ値、Ｙ値とレーザ計測によって得られたＺ値を割り当てている。
【００３７】
樹木高算出部１４は、データベース１０のＤＳＭとデータベース１１のＤＥＭとの差を求め、これをファイル１５（メモリ）のＸ，Ｙ座標に割付ける（以下樹木高グリッドデータＤｄｉという）。
【００３８】
検索範囲決定部１７は、表示部１３の画面上においてオペレータによって決められた検索範囲Ｍｉを樹冠形状指数算出部３２に設定する。
【００３９】
第１の樹冠形状指数算出部３２は、デジタル写真を基に決定された検索範囲Ｍｉに対応する樹木高グリッドデータＤｄｉをファイル１５から読み込んで、これらのグリッドデータＤｄｉにおけるアスペクト比（樹木高Ｚｃｉ/距離ｄｉ）の最大、最小値を基にして樹冠形状を明瞭にする所定の値に変更する（第１のアスペクト比更新処理３１ａ）。
【００４０】
そして、この変更された新たな最大値、新たな最小値を用いて新たな地上開度φ１´（以下第１の地上開度φ１´という）、地下開度φ２´（以下第１の地下開度φ２´という）を求め、これらの開度から算出される第１−１の樹冠形状指数φ３´を、例えば８方向で平均化して、平均化した樹冠形状指数ｅｉ（以下第１−２の樹冠形状指数ｅｉという）を得る（第１の樹冠形状指数計算処理３１ｂ）。
【００４１】
次に、局所領域Ｊｉ(例えば3×3pixel)で平滑化し、これを検索範囲Ｍｉ（例えば５ｍ）における第１−３の樹冠形状指数Ｇｅｉとしてファイル（メモリ）３７に保存する（平滑化処理３１ｃ）。
【００４２】
前述の第１−３の樹冠形状指数Ｇｅｉは、設定された樹木高グリッドデータの番号と、Ｘ、Ｙ座標と、樹木高の差Ｚと、第１−１の樹冠形状指数φ３´と、第１−２の樹冠形状指数ｅｉと、この第１−１の樹冠形状指数φ３´を得たときの地上開度φ１´と、地下開度φ２´と、検索方向毎の検索範囲以内の樹木高グリッドデータ（Ｘ、Ｙ、Ｚ）と、これらの樹木高グリッドデータまでの距離ｄｉと、設定された樹木高グリッドデータに対して高さの差Ｚ等にリンク付けされている。これらを総称して第１の樹冠形状指数情報Ｇｉという。
【００４３】
閾値算出部３５は、表示部１３に表示された画像（グリッドデータ又はオルソフォト画像又はオルソフォト画像とグリッドデータとを重ね表示）内で、予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第１−３の樹冠形状指数Ｇｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【００４４】
局所最大値フィルタ２０は、ファイル３７の第１の樹冠形状指数情報Ｇｉの各樹木高グリッドデータ（Ｘ、Ｙ、Ｚ）に対応する第１−３の樹冠形状指数Ｇｅｉに局所領域Ｆｉ（例えば３×３範囲）のフィルタを順次適用し、この結果のデータＤｆｉをファイル２７に保存する。
【００４５】
樹木頂点候補抽出部２１は、ファイル２７のデータＤｆｉとファイル３７の第１−３の樹冠形状指数Ｇｅｉとを比較し、同じ値のものを頂点候補Ｄｇｉとしてファイル２２に保存する。
【００４６】
樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の、第１−３の樹冠形状指数Ｇｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル２４に保存する。
【００４７】
第１の樹木頂点決定部２５は、ファイル２４の樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ、一致したものを仮の頂点と決定し、これをファイル２８に保存する。そして、同一樹冠部に複数の仮の頂点が存在するときは、この樹冠部で最も樹木高が高いものを単木の頂点として決定し、これをファイル２９に保存する（Ｐｉ）。
【００４８】
表示処理部１９は決定した頂点を例えば○印（△印、×印、番号等の識別）にして画面に表示する。このとき、写真画像と重ねて表示してもよい。
【００４９】
出力部３０は、ファイル２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。このとき、樹木高も表示してもよい。
【００５０】
上記のように構成された実施の形態１の樹木頂点認識装置について以下にフローチャートを用いて説明する。
【００５１】
図９は樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【００５２】
本実施の形態では、初めに初期設定を行う（Ｓ３０）。この初期設定は、例えば画面のオルソフォト画像において所定のエリアの設定、検索範囲Ｍｉ、検索方向ｋｉ等の設定である。
【００５３】
前述の検索方向ｋｉは、以下のようにして求めている。図9のS32において、任意のDdipを決定したときに、隣接するDdiを検索し、データが存在する方向の数を検索方向kiとする。例えば、任意のDdipが画像の右上角であった場合、隣接データ数は3となり、検索方向は3となる。また、任意のDdipが画像の中央であった場合、隣接データ数は8となり、検索方向は8となる。
【００５４】
本実施の形態では検索範囲Ｍｉはオルソ画像を確認し、対象林分の代表的な樹木間隔未満で設定する。
【００５５】
そして、樹冠形状指数算出部３２は、検索範囲の樹木高グリッドデータＤｄｉ（DSMとDEMデータの差分）をファイル１５から全てリードする（Ｓ３１）。
【００５６】
次に、これらの樹木高グリッドデータＤｄｉの中の１つの樹木高グリッドデータＤｄｉｐの番号を設定すると共に、検索方向内の１つの方向ｋｉｐを設定する（Ｓ３２、Ｓ３３）。本実施の形態では８方向として説明する。
【００５７】
次に、設定された番号の樹木高グリッドデータＤｄｉｐから方向ｋｉｐで検索範囲Ｍｉ内の樹木高グリッドデータＤｄｉｐ（以下グリッドデータＤｐｐｉという）を抽出する（Ｓ３４）。
【００５８】
前述の地上・地下開度を算出する検索範囲Ｍｉは対象林分に応じて変化させる必要がある。検索範囲が樹木間隔より長過ぎると隣木の頂点に対して最大傾斜角をなすため、自身が頂点として検出できない可能性がある。また、樹木間隔に対して短過ぎると樹冠内の山谷形状を検出し、一つの樹冠に対して複数の樹木頂点が抽出される場合がある。以上から検索範囲は対象林分の代表的な樹木間隔未満で設定するのが良い。
【００５９】
そして、設定された番号の樹木高グリッドデータＤｄｉｐからグリッドデータＤｐｐｉまでの間の距離ｄｉを求める（Ｓ３５）。本実施の形態では、例えば5mとして説明する。
【００６０】
次に、設定された番号の樹木高グリッドデータＤｄｉｐとグリッドデータＤｐｐｉとの垂直方向の樹木高差ΔＺｃｉ（差であることを強調するためにΔを用いている）を求める（Ｓ３６）。
【００６１】
次に、この樹木高差ΔＺｃｉを距離ｄｉで割った値（アスペクト比Ａｄｉ）を求め（Ｓ３７）、これを検索方向ｋｉｐと樹木高グリッドデータＤｄｉｐと、この番号とに対応（総称してグリッドデータ毎のアスペクト比データＶｉという）させてファイル３７に記憶する（Ｓ３８）。
【００６２】
次に、図１０に示すように、検索方向ｋｉｐが奇数か偶数かどうかを判定する（Ｓ３９ａ）。
【００６３】
ステップＳ３９ａにおいて、検索方向ｋｉｐが偶数と判定したときは、グリッドデータＤｐｐｉが７個かどうかを判定する（Ｓ３９ｂ）。ステップＳ３９ｂにおいて、グリッドデータＤｐｐｉが７個と判定したときは、処理を後述するステップＳ４１に移す。
【００６４】
また、ステップＳ３９ｂにおいて、７個ではないと判定したときは、グリッドデータＤｐｐｉを次のグリッドデータに更新して処理を図９のステップＳ３４に戻す（Ｓ３９ｃ）。
【００６５】
また、ステップＳ３９ａにおいて、検索方向ｋｉｐが奇数と判定したときは、グリッドデータＤｐｐｉが１０個かどうかを判定する（Ｓ３９ｄ）。
【００６６】
ステップＳ３９ｄにおいて、１０個と判定したときは、処理を後述するステップＳ４１に移す。
【００６７】
前述のステップＳ３９ａ、３９ｂ、３９dの判定は、例えば０．５ｍメッシュにおいて検索範囲Ｍｉが５ｍであれば、上下左右方向の場合（奇数の場合）は１０個入るので１０個かどうかを判定している。また、斜め方向の場合（偶数の場合）は７個入るので７個かどうかを判定している。
【００６８】
図１１においては、５ｍではなく１ｍの検索範囲を示して、検索範囲Ｍｉ内の検索個数を説明している。
【００６９】
また、ステップＳ４１においては、検索方向ｋｉｐが８方向になったかどうかを判定する（Ｓ４１）。
【００７０】
ステップＳ４１において、検索方向ｋｉｐが８方向全てになっていないと判定したときは検索方向を次の検索方向に更新して（Ｓ４２）、処理をステップＳ３３に戻す。
【００７１】
また、ステップＳ４１において検索方向ｋｉｐが８方向全てと判定したときは、エリア内の樹木高グリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）が最後かどうかを判定し（Ｓ４３）、最後でないときは、樹木高グリッドデータＤｄｉの番号を更新して（Ｓ４４）、処理を図９のステップＳ３２に戻す。
【００７２】
また、ステップＳ４３において、樹木高グリッドデータＤｄｉが最後と判定したときは、アスペクト比更新処理を起動する（Ｓ４５）。
【００７３】
このアスペクト比更新処理の方法の前にアスペクト比と第１の地上・地下開度（φ１´、φ２´）との関係を説明する。
【００７４】
地上・地下開度の算出には任意の点から比較する点の距離と標高差が必要となる。算出式は次式で表される。
【００７５】
Max_angle = tan-1 (最大ΔZ/d)・・・・・・式（A）
Min_angle = tan-1 (最小ΔZ/d)・・・・・・式（B）
ここで、ΔZは標高差、dは2点間の距離を表す。
【００７６】
第１の地上開度φ1´= 90−Max_angle・・・・・・式（C）
第２の地下開度φ2´= 90+ Min_angle・・・・・式（D）
方向毎に求めたアスペクト比（上記例では１０個）の中から最大及び最小のアスペクト比を抽出し、式Aと式Bから最大傾斜角と最小傾斜角を求め、式Cと式Dに当てはめて地上・地下開度（φ１´、φ２´）を算出する。但し、最大傾斜角と最小傾斜角の計算は後述するアスペクト比の更新の後となる。
【００７７】
第１の地上開度φ１´及び第２の地下開度φ２´を求めるためのアスペクト比は例えば次のように更新する（図１２参照）。
【００７８】
（1）.「ΔZ/d」が−1.7320（−６０度）以下の場合は-1.7320（-60 度）とする。
【００７９】
（樹冠端部の極大化防止）
（2）.「ΔZ/d」が-1.7320（-60 度）から0（0 度）の場合は-1.7320（-60 度）とする。
【００８０】
（樹冠山地形の強調）
（3）.「ΔZ/d」が0（0 度）から1.7320（60 度）の場合は1.7320（60 度）とする。
【００８１】
（樹冠谷地形の強調）
（4）.「ΔZ/d」が1.7320（60 度）以上の場合は「ΔZ/d」の値を据え置く。
【００８２】
（5）.「ΔZ/d」が隣接する8 方位全てがマイナスである場合「ΔZ/d」を1000倍する。
【００８３】
（樹木頂点の強調）
という処理である。
【００８４】
前述の（1）から（4）で60度を基準としてアスペクト比を更新しているが、これは現地調査でスギの樹冠形状が60度程度であったため、現地の樹冠形状に適した値として採用した。
【００８５】
すなわち、ここでは６０度と設定しているが、想定する樹冠形状によりオペレータが任意に角度を設定する。
【００８６】
更に、（５）では樹木頂点を強調するために「ΔZ/d」を1000 倍している。これにより最小角はほぼ-90度となり、最大限の樹木頂点強調が出来る。具体的には、（1）、（2）から「ΔZ/d」の変更後の値は-1.7320 となり、これを1000倍して角度を算出すると-89.9度となる。
【００８７】
このように算出された第１の地上・地下開度（φ１´、φ２´）を用いて尾根谷度の算出式（φ3 =(φ1 −φ2)/2）
と同様にして第１−１の樹冠形状指数φ3´を求める。この第１−１の樹冠形状指数φ3´を8 方位で計算し、平均化したものを第１−２の樹冠形状指数ｅｉとする。
【００８８】
以下に第１−２の樹冠形状指数ｅｉを求めるまでの処理を具体的に説明する。初めにアスペクト比更新処理について説明する。
【００８９】
アスペクト比更新処理について図１３のフローチャートを用いて説明を補充する。図１３に示すように、アスペクト比更新処理においては、ファイル３７の樹木高グリッドデータＤｄｉを設定し、検索方向ｋｉｐを設定する（Ｓ５０、Ｓ５１）。
【００９０】
次に、設定された番号の樹木高グリッドデータＤｄｉの検索方向ｋｉｐ内でのアスペクト比データViの中から値が最大なアスペクト比（ＭＡＸＡｄｉ）、値が最小のアスペクト比（ＭＩＮＡｄｉ）を抽出（読込む）する（Ｓ５２）。
【００９１】
例えば、0.5mグリッドデータを用いて、検索範囲Ｍｉが1mであるとき、ある検索方向において、アスペクト比が-2m/0.5mと-2m/1.0mであったとすると、これらの比は-4.0と-2.0になり、最大値が-2.0、最小値が-4.0となる。
【００９２】
次に、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉが上記（３）の０〜１．７３２０の範囲（角度では０度〜６０度）かどうかを判定する（Ｓ５３）。
【００９３】
ステップＳ５３において、０〜１．７３２０の範囲と判定したときは、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉを１．７３２０（ＭＡＸＡｄｉ、ＭＩＮＡｄｉ：Ａｄｂｉ）として処理をステップＳ５９に移す（Ｓ５４）。これにより、頂点から下の場合は、頂点から下の谷地形であることを強調する。
【００９４】
また、ステップＳ５３において、０〜１．７３２０の範囲（角度では０度〜６０度）ではないと判定したときは、読み込んだアスペクト比が上記（４）の１．７３２０以上かどうかを判定する（Ｓ５５）。
【００９５】
ステップＳ５５において、アスペクト比が１．７３２０以上と判定したときはそのアスペクト比（Ａｄｉ）をそのままの値（Ａｄａｉ）にして処理をステップＳ５９に移す（Ｓ５６）。
【００９６】
また、ステップＳ５５において、読み込んだアスペクト比が１．７３２０以上ではないと判定したときは、読み込んだアスペクト比が上記の（１）の−１．７３２０以下かどうかを判定する（Ｓ５７）。
【００９７】
ステップＳ５７において、−１．７３２０以下であると判定したときは、ステップＳ５９-１に処理を移す。
【００９８】
ステップＳ５９−１においては、抽出したアスペクト比を−１．７３２０（Ａｄｃｉ）にする。
【００９９】
これにより、樹冠端部と地盤との位置関係から生じるアスペクト比の極小化（樹冠形状指数では極大化）を防止する。
【０１００】
また、ステップＳ５７において、−１．７３２０以下ではないと判定したときは、抽出したアスペクト比（０〜−１．７３２０（２）(０を除く)）を−１．７３２０（Ａｄｃｉ）にする（Ｓ５９−２）。
【０１０１】
これにより、樹冠端部から上の場合は、樹冠端部から上の山地形であることを強調する。
【０１０２】
また、前述のステップＳ５４、ステップＳ５６、ステップＳ５９−２、ステップＳ５９−１のアスペクト比をファイル３７に記憶する（Ｓ５９）
次に、検索方向ｋｉｐが８かどうかを判定し（Ｓ６１）、８方向になっていないときは検索方向を更新して処理をステップＳ５１に戻す（Ｓ６２−１）。
【０１０３】
ステップＳ６２において、隣接する８方向の樹木高グリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスと判定したときは、８方向の更新後のアスペクト比を１０００倍（Ａｄｄｉ）する（Ｓ６３）。すなわち、頂点を強調する。
【０１０４】
次に、検索エリア内で樹木高グリッドデータＤｄｉが最後かどうかを判定し（Ｓ６４）、最後でないときは次の番号の樹木高グリッドデータＤｄｉに更新して処理をステップＳ５０に戻す（Ｓ６５）。
【０１０５】
すなわち、ファイル３７には、ステップＳ５９の処理によって、図１４に示すように、設定された番号の樹木高グリッドデータＤｄｉ（図１４においてはｉ＝１）と検索方向ｋｐｉと抽出したＭＡＸＡｄｉ、ＭＩＮＡｄｉと更新したアスペクト比とが対応させられて保存されている。
【０１０６】
更新したアスペクト比は、アスペクト比Ａｄｉのままの値（Ａｄａｉ）又は１．７３２０（Ａｄｂｉ）又は−１．７３２０（Ａｄｃｉ）若しくは上記3種類のいずれかの値を１０００倍したアスペクト比（Ａｄｄｉ）のいずれかである。
【０１０７】
これらを総称して本実施の形態では更新されたアスペクト比を更新アスペクト比と称する。
【０１０８】
次に、第１の樹冠形状指数算出部３２の樹冠形状指数の算出を図１６のフローチャートを用いて説明する。検索エリア内の各樹木高グリッドデータからの検索範囲において全てのアスペクト比の更新が行われると、図１６に示すように樹木高グリッドデータＤｄｉの番号を設定し（Ｓ７０）、検索方向ｋｉｐを設定する（Ｓ７１）。
【０１０９】
次に、設定されたＤｄｉ、ｋｉｐに対応する更新アスペクト比を引当、この更新アスペクト比の中から最大アスペクト比ＭＡＸＡｄｉ（アスペクト比Ａｄｉのまま又は１．７３２０又は−１．７３２０又は前記Ａｄｉを１０００倍したアスペクト比又は1732又は-1732のいずれか）と、最小アスペクト比ＭＩＮＡｄｉ（アスペクト比Ａｄｉのまま又は１．７３２０又は−１．７３２０又は前記Ａｄｉを１０００倍したアスペクト比又は1732又は-1732のいずれか）読み込む（Ｓ７２）。
【０１１０】
そして、第１の地上開度φ１´と第１の地下開度φ２´を算出する。（Ｓ７３）。
【０１１１】
第１の地上開度φ１´＝９０度−tan-1(ＭＡＸＡｄｉ)
第１の地下開度φ２´＝９０度＋tan-1(ＭＩＮＡｄｉ)
次に、この第１の地上開度φ１´と第１の地下開度φ２´とを用いて第１−１の樹冠形状指数φ３´を求める（Ｓ７４）。
【０１１２】
第１−１の樹冠形状指数φ３´＝（φ１´−φ２´）／２
次に、８方位全てに対して樹冠形状指数を求めたかどうかを判断する（Ｓ７５）。ステップＳ７５において８方位全てに対して第１−１の樹冠形状指数φ３´を求めていないと判定したときは、検索方向を更新して処理をステップＳ７１に戻す（Ｓ７６）。
【０１１３】
次に、ステップＳ７５で８方位全てに対して第１−１の樹冠形状指数φ３´を求めたと判定したときは、これを平均化（８方位）する（Ｓ７７）。
【０１１４】
そして、この平均化された第１−２の樹冠形状指数ｅｉを、樹木高グリッドデータＤｄｉを中心とした検索範囲での第１−２の樹冠形状指数ｅｉとして（Ｓ７８）、樹冠グリッドデータＤｄｉにリンク付けしてファイル３７に記憶する（Ｓ７９）。
【０１１５】
次に、検索エリア内の樹冠グリッドデータＤｄｉの番号が最後かどうかを判定し、最後でないときは樹冠グリッドデータＤｄｉの番号を更新して処理をステップＳ７０に戻す。
【０１１６】
最後に第１−２の樹冠形状指数ｅｉを局所領域Ｊｉ（例えば3×3pixel）で平滑化したものを第１−３の樹冠形状指数Ｇｅｉとして検索範囲に対応させてファイル３７に記憶する。
【０１１７】
すなわち、ファイル３７には図１４に示すように、設定された樹木高グリッドデータＤｄｉと、検索方向ｋｉｐと、更新前のＭＡＸ、ＭＩＮのアスペクト比と、更新後のＭＡＸ、ＭＩＮの更新アスペクト比とが検索方向毎に、リンクつけされて記憶されている。
【０１１８】
また、ファイル３７には図１５に示すように、設定された樹木高グリッドデータＤｄｉの番号と、そのＸ、Ｙ座標と、樹木高の差Ｚと、次の番号の樹木高グリッドデータと、このＸ、Ｙ座標と、距離ｄｉと、第１−１の樹冠形状指数φ３´と、この第１−１の樹冠形状指数φ３´を得たときの地上開度φ１´と、地下開度φ２´とが検索方向にリンクつけされ、また、８方向全てで平均化した第１−２の樹冠形状指数ｅｉ等が対応させて記憶されている。図１５には、平滑後の第１−３の樹冠形状指数Ｇｅｉを対応させて例示している。
【０１１９】
次に、頂点候補の抽出、閾値の算出方法、樹冠部の抽出処理の概略を以下に説明する。
【０１２０】
頂点候補の抽出
ここでは樹冠形状指数から樹木頂点候補を抽出する。抽出方法は次の通りである。樹冠形状指数に局所最大値フィルタを適用し、局所領域Ｆｉ内の最大値を中央画素値に置き換える。この際、局所領域Ｆｉを例えば3×3に設定した。このフィルタ処理後の指標と処理前の指標を比較し、処理前後で第１−３の樹冠形状指数Ｇｅｉが同値であったときに、その地点を「樹木頂点候補」として抽出する。
【０１２１】
局所最大値フィルタ２０は、第１−３の樹冠形状指数Ｇｅｉに対して、局所領域Fiの局所最大値フィルタを適用して、これをファイル２７に保存（Ｄｆｉ）する。
【０１２２】
次に、樹木頂点候補抽出部２１がファイル２７のデータＤｆｉを頂点候補（以下Ｄｇｉ）として抽出し、ファイル２２に保存する。具体的にはファイル３７のデータの値とファイル２７の同じ座標値のデータの値が同値の場合を「１」と、他を「０」としてファイル２２に定義する（図２０参照）。
【０１２３】
ここで、局所最大フィルタ部２０及び樹木頂点候補抽出部２１の処理について説明を補充する。
【０１２４】
局所最大フィルタ部２０（Local max filter２０ともいう）は、ファイル３７の第１の樹冠形状情報Ｇｉの第１−３の樹冠形状指数Ｇｅｉ（Ｇｅ１、Ｇｅ２・・）に局所領域Fi内で最大となる値を局所領域Fiの中央画素値に置き換える処理を順次行う。これらをＤｆｉとしてファイル２７に保存する。
【０１２５】
そして、樹木頂点候補算出部２１がデータＤｆｉとファイル３７のＧｅｉとを比較し、値が同じものを頂点候補Ｄｇｉとしてファイル２２に保存している。
【０１２６】
次に樹冠部を抽出する閾値の算出方法及び樹冠部の抽出処理を説明する。
【０１２７】
閾値以上の領域を樹冠部として抽出するが、これまで閾値は対象とする林分のオルソ画像で単木を目視確認して、同様に樹木が分割される値を閾値としていた。この手法で設定される閾値は一定値であることから樹冠部であっても閾値より小さい場合は樹冠部として抽出できなかった。この問題を解決するため、新たに画像に応じた閾値を局所領域Ｈｉ（3×3pixel や5×5pixel など）による算出で動的に決定する。
【０１２８】
初めに一つの樹冠領域を想定する。更に、この一つの樹冠領域に対して樹冠部として抽
出する領域を想定する（図１７）。ここでは尾根形状による隣木との同一化を防ぐために、樹冠領域に対して半分程度を樹木頂点が一つ含まれる樹冠部として抽出する。樹冠領域を3×3pixel から9×9pixel を想定すると樹冠部は1×1 から4×4 となり、樹冠領域に対する樹冠部の割合は平均して16％となる（表１）。
【表１】
【０１２９】
なお、0.5m グリッドデータ使用していることから、樹冠領域は1.5m×1.5m から4.5m×4.5m となり、現地で確認した樹冠領域と同等と考える。この樹冠領域を局所領域として、樹冠部の割合が16％程度になるように、局所領域の平均値（μ）に標準偏差（σ）を足した値を閾値とした。設定根拠は次の通りである。正規分布しているヒストグラムにおいて平均値の±σの範囲は全pixel の68％を占める。更に最小値から平均値の＋σの範囲を考えると84％に該当し、μ＋σ以上の範囲が16%を占める（図１８）。
【０１３０】
具体的には局所領域Ｈｉ内の平均値と標準偏差を求め、平均値と標準偏差の和を中央画素値に置き換える。第１−３の樹冠形状指数Ｇｅｉと動的閾値Ｈｇｉを比較して樹冠形状指数が大きければ、樹冠部として抽出する。
【０１３１】
次に、図１９のフローチャートを用いて具体的に樹冠部の抽出処理を説明する。初めに局所領域Ｈｉを設定する（Ｓ８０）。
【０１３２】
本実施の形態では所定エリアの樹木高グリッドデータＤｉが画像化されてオルソフォト画像に重ね表示されているとする。画面に表示された樹木高グリッドデータＤｄｉを本実施の形態では画像化されたグリッドデータＲｉと称する。
【０１３３】
この状態で、オペレータは例えば３ピクセル×３ピクセルの局所領域Ｈｉを入力する。樹冠部決定処理は、局所領域、検索エリアをメモリ（図示せず）に設定する初期設定を行う（Ｓ８０）。
【０１３４】
次に、閾値算出部が検索エリアの画像化されたグリッドデータＲｄｉを全てリードする（Ｓ８１）。
【０１３５】
次に、これらのグリッドデータＲｄｉの中からグリッドデータＲｄｎを設定し（Ｓ８２）、設定した局所領域Ｈｉに基づき、グリッドデータＲｄｎ（Ｘｄｎ、Ｙｄｎ、樹木高dn）を中心に、このグリッドデータの座標値を有する第１の樹冠形状情報Ｇｉ（３×３）をファイル３７からリードする（Ｓ８３）。
【０１３６】
このリードされた３×３の第１の樹冠形状情報Ｇｉは、第１−３の樹冠形状指数ＧｅｉをＺｉとし、グリッドデータＤｄｉのＸ軸をＤｘｉ、Ｙ軸をＤｙｉと記載する。
【０１３７】
例えば、
（Ｄｘ１、Ｄｙ１、Ｚ１）、（Ｄｘ２、Ｄｙ２、Ｚ２）、（Ｄｘ３、Ｄｙ３、Ｚ３）
（Ｄｘ４、Ｄｙ４、Ｚ４）、（Ｄｘ５、Ｄｙ５、Ｚ５）、（Ｄｘ６、Ｄｙ６、Ｚ６）
（Ｄｘ７、Ｄｙ７、Ｚ７）、（Ｄｘ８、Ｄｙ８、Ｚ８）、（Ｄｘ９、Ｄｙ９、Ｚ９）
Ｚｉ：平均化、平滑された第１−３の樹冠形状指数Ｇｅｉ
となる。
【０１３８】
そして、閾値算出部３５がこれらの第１−３の樹冠形状指数Ｇｅｉの平均値（μ）と標準偏差（σ）とを算出し（Ｓ８４）、これらの和を求め、これを設定した局所領域Ｈｉの閾値（Ｈｇｉ）として求める（Ｓ８５）。
【０１３９】
次に、閾値算出部３５は、設定したＲｄｎに対応した閾値を樹冠部抽出部２３に出力する（Ｓ８６）。樹冠部抽出部２３は第１−３の樹冠形状指数Ｇｅｉとこの閾値とを比較する（Ｓ８７）。
【０１４０】
そして、閾値（Ｈｇｉ）に対して第１−３の樹冠形状指数Ｇｅｉが大きいかどうかを判定する（Ｓ８８）。ステップＳ８８において、第１−３の樹冠形状指数Ｇｅｉの方が大きいと判定したときは、樹冠部とする（Ｓ８９）。
【０１４１】
次に、グリッドデータＲｄｉは最後かどうかを判定し（Ｓ９１）、最後でないときはグリッドデータＲｄｉを次のグリッドデータに更新し、処理をステップＳ８２に戻す（Ｓ９２）。
【０１４２】
また、ステップＳ８８において閾値（Ｈｇｉ）の方が大きいと判定したときは、樹冠部ではないと判定する（Ｓ９０）。
【０１４３】
また、樹冠部抽出部２３が閾値Ｈｇによって樹冠部のデータＤｈｉを抽出し、これをファイル２４に保存している（図２１参照）。
【０１４４】
図２１（ａ）のデータＤｈｉは、閾値を超えるＧｅｉを抽出したものであり、例えば番号１のエリアの第１−３の樹冠形状指数Ｇｅｉは、（Ｘｈ５、Ｙｈ５）を「１」、（Ｘｈ６、Ｙｈ６）を「１」・・（Ｘｈ１０、Ｙｈ１０）を「１」として示している。
【０１４５】
図２１の（ｂ）は図２１（ａ）を具体的に示すものであり、閾値を超えた画素は「１」とする。「１」が連続している範囲が樹冠部である。
【０１４６】
そして、ファイル２４の中の頂点候補ＤｈｉのＸｈｉ、Ｙｈｉをポリゴン化（座標Ｘｈｉ、Ｙｈｉと色情報）して表示処理部１９によってラベルリング（纏まったエリアを（樹冠部）にラベルを振る）させる（図２６参照）
次に、樹木頂点決定部２５は、ファイル２４のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・）とファイル２２の頂点候補Ｄｇｉ（Ｄｇ１、Ｄｇ２・・）とを重ね合わせ一致したものを仮の頂点と決定する。そして、樹冠部のエリア内に複数の仮の頂点が存在するときは、この樹冠部のエリアで最も樹木高が高いものを単木の頂点として決定し、これをファイル２９に保存する（Ｐｉ）。表示処理部１９は決定した頂点を画面に表示する（図２３参照）
前述の頂点の決定について説明を補充する。例えば、図２２の（ａ）に示すように、ファイル２４に第１列目が「００１１００１１１００」、第２列目が「０１１１１０１１１１０」、第３列目が「００１１０００１０００」、第４列目が「０００００１０００１１」という具合に頂点上部域が得られた場合に、図２２（ｂ）に示すように、「１」が隣接するまとまり（「０」が区切り）を１つの樹冠部としてまとまり毎に通し番号を割り振っていく。
【０１４７】
そして、例えば、同一高の頂点が樹冠内に複数の頂点が存在することになると、一度、ファイル２８のデータＤＱｉを表示させて、オペレータによって指定させ、指定されたものを頂点Ｐｉと決定してファイル２９に保存する。
【０１４８】
このような処理によって第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出結果は、図８、図２３に示すように頂点が求められている。
【０１４９】
また、表２に尾根谷度と第１−３の樹冠形状指数Ｇｅｉの樹木頂点抽出結果の比較を示す。
【表２】
【０１５０】
表２から、抽出精度は同程度であるが、抽出率及び正解率が向上していることが分かる。以上から第１−３の樹冠形状指数Ｇｅｉの有意性が確認できた。
【０１５１】
また、図２４の尾根谷度画像では林道に接した樹木は極大化して画像が明るくなっているが図２５の樹冠形状指数画像では同箇所が暗くなり、極大化を解消していることがわかる。
【０１５２】
次に、尾根谷度及び第１−３の樹冠形状指数による樹木頂点の抽出数、判読による樹木頂点数と現地で樹冠上層部に到達していると確認した樹木本数を表３にまとめる。
【表３】
【０１５３】
表から、尾根谷度による樹木頂点抽出はP8-2 を除いて3〜4 割の程度の抽出率であったのに対して、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出は7〜9 割の抽出率であった。また、空中写真からの情報では密度の高い森林において樹木の抽出に限界があることがわかっている。このことから、同様に上空からの情報であるレーザデータを用いた樹木頂点抽出では判読による樹木頂点数を検証値として考えることもできる。
【０１５４】
そこで、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出と判読による樹木頂点の割合を算出した。その結果、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出数は判読の樹木頂点数の8 割強〜10 割強であった。以上から、第１−３の樹冠形状指数Ｇｅｉによる樹木頂点抽出は尾根谷度による抽出よりも倍近い抽出率が得られ、なおかつ、目視判読による樹木頂点抽出とほぼ同レベルで行えることが検証できた。
【０１５５】
また、本実施の形態は、図２７に示すようなコンピュータシステムVPSによって実現するのが好ましい。ＶＰＳ１は、ワークステーション、プロセッサ、マイクロコンピュータ、ロジック、レジスタ等の適宜な組み合わせからなる中央情報処理装置（CPU）と、この中央情報処理装置に必要な制御・操作情報を入力するキーボード（KB）、マウス、対話型ソフトスイッチ、外部通信チャネル等を含む情報入力部と、中央情報処理装置から出力された情報を広義な意味で表示・伝送するディスプレイ、プリンタ、外部通信チャネル等を含む情報出力部と、中央情報処理装置に読み込まれるオペレーティングシステム、アプリケーションプログラム等の情報が格納されたロム（ＲＯＭ）と、中央情報処理装置で随時処理すべき情報及び中央情報処理装置から随時書き込まれる情報を格納するラム（ＲＡＭ）等とを備える。ＲＯＭ、ＲＡＭを適宜統合、細分化することは差し支えない。
【０１５６】
前述のアプリケーションプログラムは、本実施の形態の、樹木高算出部１４と、検索範囲決定部１７と、局所領域設定部３６と、第１の樹冠形状指数算出部３２と、閾値算出部３５と、局所最大フィルタ２０（Ｌｏｃａｌ max ｆｉｌｔｅｒともいう）と、樹木頂点候補抽出部２１と、樹冠部抽出部２３と、第１の樹木頂点決定部２５と、出力部３０等である。
【０１５７】
なお、上記実施の形態では、画面において決定した検索エリアにおける樹木の各樹木頂点を求めるとしたが、検索エリアではなく、樹木画像上で軌跡Ａｉを描き、この軌跡Ａｉ上における樹木の頂点を抽出してもよい。
【０１５８】
つまり、軌跡Ａｉの樹木高グリッドデータを抽出し、この樹木高グリッドデータに対応する第１−３の樹木形状指数Ｇｅｉを読み出し、この樹木形状指数Ｇｅｉによって頂点を決定する。
【０１５９】
＜実施の形態２＞
スギ、ヒノキ人工林は通常3000 本/ha 植栽され、林齢17 年程度で第1 回間伐が行われる。伐採率は3 割程度であるが、自然枯死を併せると間伐後は1950 本/ha 程度となる。
【０１６０】
その後、6〜8 年毎に第2 回、第3 回の間伐（間伐率は共に約3 割）を行い、950 本/ha 程度の立木密度にして収穫を行う。現在日本のスギ、ヒノキの植林は少なく、そのほとんどは林齢20 年以上であり、第1 回の間伐は実施されていることと想定される。
【０１６１】
このことから、第2 回の間伐が終了しているか判断するためにも2000 本/ha 程度の立木密度の森林を把握することが求められる。
【０１６２】
本実施の形態２は2000 本/ha 程度の立木密度の密な森林（例えばスギ、ヒノキ等）においても樹木頂点の抽出精度を向上できる。
【０１６３】
本実施の形態２では、抽出精度の向上のために、「地上開度、地下開度の算出の見直し」と「局所加重平均」と、「樹木頂点が隣接する場合に所定の樹木頂点を除去する処理」、「森林領域の抽出」等を備える。
【０１６４】
（地上開度φ１´、地下開度φ２´の算出の見直し）
実施の形態１は疑似樹冠モデルという考えをもとに、DHM を用いて算出した地上・地下開度を樹冠モデルに合わせて数値を置き換えていた。このとき、樹冠端部から樹木頂点の角度を60°と設定していたが、本実施の形態２では現地で確認したスギ・ヒノキの樹形に合うように70°へと変更する。
【０１６５】
また、実施の形態１では林縁部の極大化を防ぐために、林縁部から地表を見下ろしたときに−60°以上の角度となる場合、林縁部の値を−60°に置き換えていた。しかし、依然として林縁部は凸の形状となり、第１−３の樹冠形状指数Ｇｅｉは大きく、一方、隣接する樹冠の境界に位置する場合は、凹の形状となり、第１−３の樹冠形状指数Ｇｅｉは小さくなる傾向にある。
【０１６６】
また、どちらも樹冠端部に位置するという点では一致しているものの、第１−３の樹冠形状指数Ｇｅｉに大きな違いが生じる。そこで、この違いを解消するため、新たな樹冠形状の−70°以上の角度となるときに林縁部と規定し、その値を70°(隣接する樹冠の境界に位置する場合の角度)に置き換える。なお、ここではDHM が3m 以下の非森林領域は樹冠形状指数の最小値（-70）とした。詳細については後述する。
【０１６７】
（局所加重平均）
実施の形態１では、算出された地上・地下開度（φ１、φ２）を用いて尾根谷度の算出式（φ3 = (φ1 −φ2 )/ 2）と同様な計算式で第１−１の樹冠形状指数φ3´ を求め、この第１−１の樹冠形状指数φ3´ を8 方位で計算して、平均化した第１−２の樹冠形状指数ｅｉを更に局所領域で平滑化処理して第１−３の樹冠形状指数Ｇｅｉとしていた。疎林において、この平滑化は複雑な樹冠形状を滑らかにして、同一樹冠内に生じる複数の樹木頂点を絞る効果があった。
【０１６８】
しかしながら、密林においては上記作用により樹木頂点の抽出数を減らすという結果になり、抽出精度が低くなる場合がある。隣接する樹木頂点が平滑化により樹木頂点数が減少する例を図２８に示す。図の上部分は樹冠形状指数の画像を数値化したものであり、下部分は横線部（２８ａ、２８ｂ）の断面形状を示している。
【０１６９】
図２８から、同一局所領域内に樹木頂点が複数含まれる場合、[１]平滑化により樹木頂点が減少する、[２]樹木頂点の位置座標が変わることが問題であるとわかる。このため、本実施の形態２では後述する局所加重平均によってこの問題を解決する。
【０１７０】
（森林領域の抽出）
樹冠形状指数は地表面（DSM）と地盤面（DEM）の差分（DHM）をもとにわずかな凹凸からも樹冠モデルを仮想的に作成するものである。一方で、樹木が存在しないギャップや林道などの領域においても樹冠モデルを作成してしまい、樹木が存在しない領域でも樹木頂点の抽出がなされてしまう。
【０１７１】
そこで、樹木が存在しない領域において樹木頂点の抽出を行わないように、森林領域の抽出を行った（図２９参照）。森林領域の抽出にはDHM を用い、下層植生を含まないようにここでは一例としてDHMが3m以上の領域を森林領域とした。
【０１７２】
図３０は本実施の形態２の樹木頂点認識装置の概略構成図である。図３１は本実施の形態２の樹木頂点認識装置を説明する概略フローチャートである。初めにこの概略フローチャートを説明する。
【０１７３】
図３１に示すように、航空レーザ計測により取得されたデータ（ＤＳＭ：Digital Surface Model）を読み込む（Ｓ３０）。また、地盤のグリッドデータ（ＤＥＭ：Digital Elevation Model）を読み込む（Ｓ３１）。
【０１７４】
そして、DSMとDEMデータの差分Ｚｃｉ（Ｚｃｉ、Ｚｃｉ・・）を樹木高グリッドデータＤｄｉとして求めて（Ｓ３２）、樹冠形状の再現処理を行う（Ｓ３３）。
【０１７５】
この樹冠形状の再現処理は、樹木高のグリッドデータＤｄｉ（Ｄｄ１、Ｄｄ２、・・・）に対して樹木間隔程度（例えば１．５ｍ〜５ｍ範囲（図７を参照））を設定する（Ｓ３３ａ）。
【０１７６】
次に、各グリッドデータと自身のデータとのアスペクト比を求め、このアスペクト比を一定条件下で新たなアスペクト比（スギ、ヒノキの樹冠形状に合わせて７０°に変更）に合わせて変更し、変更したアスペクト比（（ΔＺｉ／ｄｉ）＝Ａｄｉ）を用いて第２の地上開度φ１´及び第２の地下開度φ２´を算出する（Ｓ３３ｂ）。
【０１７７】
本実施の形態２では現地で確認したスギ・ヒノキの樹形に合うように70°へと変更する。また、実施の形態１では、林縁部の極大化を防ぐために、林縁部から地表を見下ろしたときに−60°以上の角度となる場合、林縁部の値を−60°に置き換えていた。
【０１７８】
しかし、依然として林縁部は凸の形状となり、樹冠形状指数は大きく、一方、隣接する樹冠の境界に位置する場合は、凹の形状となり、樹冠形状指数は小さくなっていた。
【０１７９】
どちらも樹冠端部に位置するという点では一致しているものの、樹冠形状指数に大きな違いがあった。そこで、この違いを解消するため、新たな樹冠形状の−70°以上の角度となるときに林縁部と規定し、その値を70°(隣接する樹冠の境界に位置する場合の角度)に置き換える。
【０１８０】
そして、第２の地上開度φ１´及び第２の地下開度φ２´を尾根谷度の算出式により、第２−１の樹冠形状指数φ３´を算出してこれを８方位で平均化した第２−２の樹冠形状指数ｅｉを得た後で、ガウシアンフィルタを用いた平滑化（図３２参照）を行った第２−３の樹冠形状指数ｅｅｉを得る。次に、森林領域（例えばDHMが3m以上の領域）を規定し、この森林領域における第2-3の樹冠形状指数eeiはそのままの値とし、非森林領域（例えばDHMが3m未満の領域）の樹冠形状指数を任意の数値（例えば「-70」）とする第２−４の樹冠形状指数Ｅｅｉを求め（Ｓ３３ｃ）（図２９右図参照）、樹冠形状指数画像を得る（Ｓ３３ｄ）。
【０１８１】
この第２−４の樹冠形状指数Ｅｅｉ（単に樹冠形状指数Ｅｅｉという場合もある）と実施の形態１の第１−３の樹冠形状指数Ｇｅｉとを総称して樹冠形状指数ＥＧｉという。
【０１８２】
前述の開度の算出で、樹木間隔を樹木間隔程度としたが、これは検索範囲が長すぎると過小抽出となり、短すぎると過大抽出となるためである。一例として、表４に示すように立木密度から樹木間隔を求めるのが好ましい。
【表４】
【０１８３】
但し、Ｈ１〜Ｈ５は検証ポイント（所定の領域を有する）を示す。
【０１８４】
次に、樹木頂点の算出処理を行う（Ｓ３４）。
【０１８５】
この樹木頂点の算出処理は、第２−４の樹冠形状指数Ｅｅｉと動的閾値との比較により樹冠部のデータＤｈｉを抽出して（Ｓ３４ａ）、樹冠部画像を得る（Ｓ３４ｂ）。この樹冠部域の抽出を行う動的閾値の局所領域サイズは樹木間隔程度である。
【０１８６】
局所領域を樹木間隔程度とするのは、樹木毎の第２−４の樹冠形状指数Ｅｅｉに応じて樹冠の上部領域を抽出するためである。
【０１８７】
動的閾値は局所領域の平均値（μ）に標準偏差（σ）を足した値であり、この閾値に
より1 本の樹木の樹冠領域に対して約16％の中央部（樹冠の上部領域）を抽出することができる。
【０１８８】
また、樹冠部域の抽出は上記の通り局所領域で閾値を求め、この値と局所領域の中心画素の第２−４の樹冠形状指数Ｅｅｉとを比較して第２−４の樹冠形状指数Ｅｅｉが大きければ樹冠部域のデータＤｈｉ（ＤＱｉ）として抽出する。
【０１８９】
表9に樹木間隔と各プロットの局所領域サイズを示す。
【表９】
この局所領域サイズを用いて抽出した樹冠部域を図３４に示す。
【０１９０】
そして、樹冠部のデータのデータＤＱｉに対応する第２−４の樹冠形状指数Ｅｅｉの中から樹木冠形状指数が最高値となる第２−４の樹冠形状指数Ｅｅｉを検索し、これを樹木頂点Pipとする（Ｓ３４ｃ）。
次に、検索した樹木頂点の隣接程度を判断し、隣接していると判断したときは、第２−４樹冠形状指数Ｅｅｉの大小関係を基に隣接する樹木頂点のデータＰｉｐの第２−４の樹冠形状指数Ｅｅｉの大小関係をもとに、樹木頂点のデータＰｉｐを除去して（Ｓ３４ｄ）、樹木頂点Ｐｉとした画像（図３５参照）を得る（Ｓ３４ｅ）。
【０１９１】
次に、図３０の構成を説明する。
【０１９２】
（図３０の説明）
図３０においては図１と同一の符号のものについては説明を省略する。図３０に示すように、本実施の形態２の樹木頂点認識装置は、第２の樹冠形状指数算出部４１と、第２の樹木頂点決定部４３（隣接樹木除去処理を含む）等を備える。
【０１９３】
第２の樹冠形状指数算出部４１は、実施の形態１と同様な樹冠形状指数の算出処理（図示せず）を備えると共に、新たな第２のアスペクト比更新処理４１ａと、第２の樹冠形状指数計算処理４１ｂと、ガウシアンフィルタ処理４１ｃと、森林領域抽出処理４１ｄ等を備えている。
【０１９４】
第２のアスペクト比更新処理４１ａは、樹木高グリッドデータＤｄｉ（Ｘ、Ｙ、ＤＳＭとＤＥＭの差（ＤＳＭ））をファイル１５から全て読み込んで、これらのグリッドデータＤｄｉにおけるアスペクト比（樹木高Ｚｃｉ/距離ｄｉ）を求める。このアスペクト比は、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９５】
そして、図３７、図３８に示す様にアスペクト比を更新する。このアスペクト比が−70°以上の角度となるときに、その値を70°に置き換える。また、−７０°から０度の場合は−７０°に置き換える。この新たな更新アスペクト比は、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。また、本実施の形態２の新たなアスペクト比の更新については詳細に後述する。
【０１９６】
第２の樹冠形状指数計算処理４１ｂは、アスペクト比更新処理が終わると、３×３の領域毎にこの変更された新たな最大値、新たな最小値を用いて新たな第２の地上開度φ１´、第２の地下開度φ２´を求め、これらの開度から第２−１の樹冠形状指数φ３´を求めて、８方向で平均化して、第２−２の樹冠形状指数ｅｉを順次求める。
【０１９７】
ガウシアンフィルタ処理４１ｃは、第２−２の樹冠形状指数ｅｉの算出処理が終わると、３×３の領域毎の第２−２の樹冠形状指数ｅｉに対して局所加重平均化処理を行った第２−３の樹冠形状指数ｅｅｉを得る。この第２−３の樹冠形状指数ｅｅｉはグリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９８】
森林領域抽出処理４１ｄは、ガウシアンフィルタ処理後の第２−３の樹冠形状指数ｅｅｉの関連する樹木高グリッドデータのDHMが例えば３ｍ以上の場合に、森林領域の第２−４の樹冠形状指数Ｅｅｉとして抽出する。この第２−４の樹冠形状指数Ｅｅｉは、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０１９９】
閾値算出部３５は、表示部１３に表示された画像（グリッドデータ又はオルソフォト画像又はオルソフォト画像とグリッドデータとを重ね表示）内で、予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第２−４の樹冠形状指数Ｅｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【０２００】
樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の第２−４の樹冠形状指数Ｅｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・、図２２（ａ）参照）としてファイル２４に保存する。
【０２０１】
第２の樹木頂点決定部４３は、メモリ２４のデータＤｈｉの隣接するまとまり毎に番号を割り振ったデータＤＱｉ（樹冠形状指数のデータ）をメモリ２８に生成し、まとまり毎の番号の内で最大の第２−４の樹冠形状指数Ｅｅｉを持つ点を樹木頂点Pipとする。
【０２０２】
さらに、樹木頂点Pipに対して例えば3×3の局所領域を設定し、順次この局所領域を移動させ、局所領域の中央に樹木頂点Pipが位置するときに局所領域内に他の樹木頂点Pipが存在するかどうかを判定する。存在している場合は、第２−４の樹冠形状指数Ｅｅｉの大小を比較し、中央の樹木頂点Pipの樹冠形状指数Ｅｅｉが小さい場合、この樹木頂点を消去する。
【０２０３】
そして、隣接樹木消去処理が終わったデータを樹木頂点Ｐｉとしてメモリ２９に記憶する。
【０２０４】
出力部３０は、ファイル２９の頂点Ｐｉの座標（ＸＹ）に樹木番号を付けてプリンタ（図示せず）又は画面に表示する。
【０２０５】
上記のように構成された実施の形態の樹木頂点認識装置について以下にフローチャートを用いて説明する。
【０２０６】
本実施の形態２では、上記の実施の形態１と同様な樹冠立体角算出処理によって、既に検索方向ｋｉｐ毎のグリッドデータＤｐｐｉと樹木高グリッドデータＤｄｉｐとの更新アスペクト比Ａｄｉが求められてファイル３７に記憶され、アスペクト比更新処理が起動したとする。
【０２０７】
図３６は実施の形態２のアスペクト比更新処理を説明するフローチャートである。
【０２０８】
実施の形態２のアスペクト比の更新は図３７に示すように更新する。
【０２０９】
つまり、
（6）.「ΔZ/d」が-2.7474（-70 度）以下の場合は2.7474（70 度）とする。
【０２１０】
（樹冠端部の極大化防止）
（7）.「ΔZ/d」が-2.7474（-70 度）から0（0 度）の場合は-2.7474（-70 度）とする。
【０２１１】
（樹冠山地形の強調）
（8）.「ΔZ/d」が0（0 度）から2.7474（70 度）の場合は2.7474（70 度）とする。
【０２１２】
（樹冠谷地形の強調）
（9）.「ΔZ/d」が2.7474（70 度）以上の場合は「ΔZ/d」の値を据え置く。
【０２１３】
（樹木頂点の強調）
（10）.「ΔZ/d」が隣接する8 方位全てがマイナスである場合「ΔZ/d」を1000倍する。
【０２１４】
という具合に更新する。
【０２１５】
アスペクト比更新処理は、図３９に示すようにファイル３７の樹木高グリッドデータＤｄｉを設定し、検索方向ｋｉｐを設定する（Ｓ１００、Ｓ１０１）。
【０２１６】
次に、設定された番号の樹木高グリッドデータＤｄｉの検索方向ｋｉｐ内でのアスペクト比データViの中から値が最大なアスペクト比（ＭＡＸＡｄｉ）、値が最小のアスペクト比（ＭＩＮＡｄｉ）を抽出する（Ｓ１０２）。
【０２１７】
次に、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉが上記（８）の０〜２．７４７４の範囲（角度では０度〜７０度）かどうかを判定する（Ｓ１０３）。
【０２１８】
ステップＳ１０３において、０〜２．７４７４の範囲と判定したときは、読み込んだアスペクト比ＭＡＸＡｄｉ、ＭＩＮＡｄｉを２．７４７４（ＭＡＸＡｄｉ、ＭＩＮＡｄｉ：Ａｄｂｉ）として処理をステップＳ１１１に移す（Ｓ１０４）。これにより、頂点から下の場合は、頂点から下の谷地形であることを強調する。
【０２１９】
また、ステップＳ１０３において、０〜２．７４７４の範囲（角度では０度〜７０度）ではないと判定したときは、読み込んだアスペクト比が上記（９）の２．７４７４以上かどうかを判定する（Ｓ１０５）。
【０２２０】
ステップＳ１０５において、アスペクト比が２．７４７４以上と判定したときはそのアスペクト比（Ａｄｉ）をそのままの値（Ａｄａｉ）にして処理をステップＳ１１１に移す（Ｓ１０６）。
【０２２１】
また、ステップＳ１０５において、読み込んだアスペクト比が２．７４７４以上ではないと判定したときは、読み込んだアスペクト比が上記の（６）の−２．７４７４以下かどうかを判定する（Ｓ１０７）。
【０２２２】
ステップＳ１０７において、−２．７４７４以下であると判定したときは、ステップＳ１０９に処理を移す。ステップＳ１０９においては、抽出したアスペクト比を＋２．７４７４（Ａｄbｉ）にする。
【０２２３】
これにより、樹冠端部と地盤との位置関係から生じるアスペクト比の極小化（樹冠形状指数では極大化）を防止する。
【０２２４】
また、ステップＳ１０７において、−２．７４７４以下ではないと判定したときは、抽出したアスペクト比（０〜−２．７４７４（７）(０を除く)）を−２．７４７４（Ａｄｃｉ）にする。これにより、樹冠端部から上の場合は、樹冠端部から上の山地形であることを強調する。
【０２２５】
また、前述のステップＳ１０４、ステップＳ１０６、ステップＳ１１０、ステップＳ１０９のアスペクト比をファイル３７に記憶する（Ｓ１１１）
次に、検索方向ｋｉｐが８かどうかを判定し（Ｓ１１２）、８方向になっていないときは検索方向を更新して処理をステップＳ１０１に戻す（Ｓ１１３）。
【０２２６】
ステップＳ１１３において、隣接する８方向の樹木高グリッドデータのアスペクト比が、水平軸を基準として全てマイナスと判定したときは、８方向の更新後のアスペクト比を１０００倍（Ａｄｄｉ）する（Ｓ１１５）。すなわち、頂点を強調する。
【０２２７】
次に、検索エリア内で樹木高グリッドデータＤｄｉが最後かどうかを判定し（Ｓ１１６）、最後でないときは次の番号の樹木高グリッドデータＤｄｉに更新して処理をステップＳ１００に戻す（Ｓ１１７）。
【０２２８】
すなわち、ファイル３７には、ステップＳ１０９、Ｓ１１０の処理によって、図３８に示すように、設定された番号の樹木高グリッドデータＤｄｉと検索方向ｋｐｉと抽出したＭＡＸＡｄｉ、ＭＩＮＡｄｉと更新アスペクト比とが対応させられて保存されている。
【０２２９】
更新アスペクト比は、アスペクト比Ａｄｉのままの値（Ａｄａｉ）又は２．７４７４（Ａｄｂｉ）又は−２．７４７４（Ａｄｃｉ）若しくは−２．７４７４を１０００倍したアスペクト比（Ａｄｄｉ）のいずれかである。
【０２３０】
従って、第２の樹冠形状指数計算処理４１ｂは、Ｄｄｉ、ｋｉｐに対応する更新アスペクト比を引当、この更新アスペクト比の中から最大アスペクト比ＭＡＸＡｄｉ（アスペクト比Ａｄｉのまま又は２．７４７４又は−２．７４７４又は−２７４７．４のいずれか）と、最小アスペクト比ＭＩＮＡｄｉ（アスペクト比Ａｄｉのまま又は２．７４７４又は−２．７４７４又は−２７４７．４のいずれか）読み込むことになる。
【０２３１】
（第２の樹冠形状指数計算処理）
そして、第２の樹冠形状指数計算処理４１ｂは、第２の地上開度φ１´と第２の地下開度φ２´を算出する。
【０２３２】
第２の地上開度φ１´＝９０度−tan-1(ＭＡＸＡｄｉ)
第２の地下開度φ２´＝９０度＋tan-1(ＭＩＮＡｄｉ)
次に、この第２の地上開度φ１´と第２の地下開度φ２´とを用いて第２−１の樹冠形状指数φ３´を求める。
【０２３３】
第２−１の樹冠形状指数φ３´＝（φ１´−φ２´）／２
次に、８方位全てに対して、これを平均化（８方位）し、この平均化された第２−１の樹冠形状指数φ３´を、樹木高グリッドデータＤｄｉを中心とした検索範囲での第２−２の樹冠形状指数ｅｉとして、樹冠グリッドデータＤｄｉにリンク付けしてファイル３７に記憶する。
【０２３４】
最後に第２−２の樹冠形状指ｅｉを検索範囲に対応させてファイル３７に記憶する。
【０２３５】
次に、本実施の形態２で用いるガウシアンフィルタ処理４１ｃについて説明する。
【０２３６】
この手法は局所平均（これまでの平滑化）と異なり、対象局所領域に重みを持たせ、画像をより鮮明にすることができる。局所平均と局所加重平均の例を図４０に示す。図４０（ａ）が局所平均であり、図４０（ｂ）が局所加重平均の処理結果である。
【０２３７】
つまり、本実施の形態ではガウシアンフィルタという局所加重平均を用いて第２−２の樹冠形状指数ｅｉを平滑化する。
【０２３８】
ガウシアンフィルタは下記の式（E）のガウス関数で加重を決定する。
【数１】
【０２３９】
ここで、G は各画素の加重、i,j は縦横それぞれの中心画素からの距離（画素）、νは変数を表す。
【０２４０】
ガウシアンフィルタは変数νを任意で変えることで、画像の平滑化の度合いを変化させることができるという利点を持っている。
【０２４１】
一般的に局所領域の大きさが(2N＋1)×(2N＋1)
の場合、νはN/2 とするのが目安である。また、νの値が小さいとより原画像に近く、νの値が大きいと局所平均画像に近づく特徴がある。図４４にガウシアンフィルタを用いた平滑化の例を示す。なお、νは目安となっているN/2 とした。図から上記で述べた局所平均による平滑化の2つの問題点を解決できていることから、本実施の形態２ではガウシアンフィルタを用いた平滑化を採用する。νは、０．５〜０．７（図３２参照）とする（好ましくは０．７）。
【０２４２】
ガウシアンフィルタ処理４１は、このような第２−３の樹冠形状指数ｅｅｉを、メモリ３７に記憶し、森林領域抽出処理４１ｄが以下の処理を行う。
【０２４３】
（森林領域抽出処理）
次に、森林領域抽出処理を説明する。
【０２４４】
樹冠形状指数は地表面（DSM）と地盤面（DEM）の差分（DHM）をもとにわずかな凹凸からも樹冠モデルを仮想的に作成するものである。一方で、樹木が存在しないギャップや林道などの領域においても樹冠モデルを作成してしまい、樹木が存在しない領域でも樹木頂点の抽出がなされてしまう。
【０２４５】
そこで、樹木が存在しない領域において樹木頂点の抽出を行わないように、森林領域の抽出を行った（図２９参照）。森林領域の抽出にはDHM を用い、下層植生を含まないようにここでは一例としてDHMが3m以上の領域を森林領域とした。そして、森林領域の第2-3の樹冠形状指数eeiはそのままの値を用い、非森林領域（例えばDHMが3m未満の領域）の樹冠形状指数を任意の数値（例えば「-70」）にしたものを第２−４の樹冠形状指数Ｅｅｉとした。この第２−４の樹冠形状指数Ｅｅｉは、グリッドデータＤｄｉに関連付けられてメモリ３７に記憶される。
【０２４６】
このため、ファイル３７には図３９に示すように、設定された樹木高グリッドデータＤｄｉと、検索方向ｋｉｐと、更新前のＭＡＸ、ＭＩＮのアスペクト比と、更新後のＭＡＸ、ＭＩＮのアスペクト比とが検索方向毎に、リンクつけされて記憶されている。
【０２４７】
また、ファイル３７には図３９に示すように、設定された樹木高グリッドデータＤｄｉの番号と、そのＸ、Ｙ座標と、樹木高の差Ｚと、次の番号の樹木高グリッドデータと、このＸ、Ｙ座標と、距離ｄｉと、第２―１の樹冠形状指数φ３´と、この第２―１の樹冠形状指数φ３´を得たときの第２の地上開度φ１´と、第２の地下開度φ２´とが検索方向にリンクつけされ、また、これらの８方向全ての平均化された第２−２の樹冠形状指数ｅｉ及びガウシアンフィルタ後の第２−３の樹冠形状指数ｅｅｉ、森林領域抽出後の第２−４の樹冠形状指数Ｅｅｉ等が対応させて記憶されている。これらを総称して第２の樹冠形状指数情報Ｇｉｉともいう。
【０２４８】
次に、閾値算出部３５が予め設定された局所領域Ｈｉ（例えば３×３ピクセル）の第２−４の樹冠形状指数Ｅｅｉの平均値μと標準偏差σとの和を閾値Ｈｇｉとして求め、これを樹冠部抽出部２３に渡す。
【０２４９】
次に、樹冠部抽出部２３は、閾値算出部３５で求められた閾値Ｈｇｉより大きい値の第２−４の樹冠形状指数Ｅｅｉを抽出し、これを樹冠部のデータＤｈｉ（Ｄｈ１、Ｄｈ２・・・）としてファイル２４に保存する。
【０２５０】
次に、第２の樹木頂点決定部４３は、メモリ２４のデータＤｈｉの隣接するまとまり毎に番号を割り振ったデータＤＱｉ（樹幹形状指数のデータ）をメモリ２８に生成する（図４６（ａ））。
【０２５１】
このデータＤＱｉに対して以下の樹木頂点抽出処理を行う。
【０２５２】
データＤＱｉは樹冠領域毎（まとまり毎）に重複しない番号が割り振られる。一つの番号につき、複数のデータDdiが存在する場合、この番号の内で最大の樹冠形状指数を持つ点を樹木頂点Pipとする。
【０２５３】
さらに、樹木頂点が隣接し、過剰抽出となる場合があるため、隣接樹木頂点の除去処理を行う。
【０２５４】
樹木頂点抽出結果（上記手法により抽出された樹木頂点Pipの抽出精度）からスギ林だけでなく、ヒノキ林においても精度良く抽出する手法が開発できたと考える。なお、樹冠部域画像と樹木頂点抽出結果を参照すると、抽出樹木頂点が隣接している場合に若干多めに樹木頂点を抽出する傾向にあるが、この点に対しては以下の手法により解決できる。
【０２５５】
図３５に1 つの樹冠に複数の樹木頂点が存在している状況を示す。これは、樹冠部域が1 つの樹冠に対して1 つの樹冠部域が抽出されず、樹冠部が分割されたために生じたと考えられる0.5m グリッドにおいて樹木頂点が斜めに隣接する場合、2 点間の距離は0.7m であり、立木密度を考慮するとその間隔は短い。このことからも1 つの樹冠から過剰に樹木頂点が抽出されていることがわかる。
【０２５６】
この問題を解決するために、樹木頂点の抽出でも使用した樹冠形状指数の大小関係を利用して、隣接する樹木頂点を除去した。具体的には、Pipに対して例えば3×3の局所領域を設定し、順次この局所領域を移動させ、局所領域の中央に樹木頂点Pipが位置するときに局所領域内に他の樹木頂点Pipが存在するかどうかを判定し、存在している場合は、樹冠形状指数Ｅｅｉの大小を比較し、中央の樹木頂点Pipの樹冠形状指数Ｅｅｉが小さい場合、この樹木頂点を消去する。
【０２５７】
図４５は隣接樹木除去処理を説明するフローチャートである。
【０２５８】
隣接樹木処理は、樹木頂点のデータPipに対して３×３のフィルタをかけて、このフィルタの中央の位置ＪＤｉを設定する（Ｓ２０１、図４６（ｂ））。
【０２５９】
次に、中央の位置ＪＤｉの隣接8方向の領域ＪＤｒｉに樹木頂点のデータPipがあるかどうかを判定する（Ｓ２０２）。
【０２６０】
ステップＳ２０２において、データＰｉｐがあると判定したときは、隣接8方向の領域ＪＤｒｉにデータＰｉｐに対応する第２−４の樹冠形状指数Ｅｅｉと中央の位置ＪＤｉの第２−４の樹冠形状指数Ｅｅｉとを比較して値の大小を判定する（Ｓ２０３、Ｓ２０４）。
【０２６１】
ステップＳ２０４において、ＪＤｉの樹冠形状指数Ｅｅｉが小さいと判定したときは、ＪＤｉの樹木頂点のデータPipを「0」とする（Ｓ２０６）。
【０２６２】
また、ステップＳ２０４において、ＪＤｉの樹冠形状指数Ｅｅｉが大きいと判定したときは、樹木頂点のデータPipがメモリ２８に他にあるかどうかを判定する（Ｓ２０７）。他にあると判定したときは、ＪＤｉを次の番号に更新して処理をステップＳ２０１に移す（Ｓ２０９）。
【０２６３】
図４６（ｂ）においては、３×３フィルタの中央に「２」が存在し、斜め方向に「１」が存在している。「２」に対応する第２−４の樹冠形状指数が、「１」に対応する第２−４の樹冠形状指数より小さい場合は、図４６（ｃ）に示すように「２」は消去される。
【０２６４】
また、ステップＳ２０７において、データＰｉｐが他にないと判定したときは、メモリ２８の全ての樹木頂点のデータPipを樹木頂点Ｐｉとしてメモリ２９に記憶する（Ｓ２０８）。
【０２６５】
次に、樹木頂点消去について説明を補充する。
図４４ に樹木頂点を除去した結果を示す。図に示す通り、同一樹冠内の樹木頂点が除去され、1つの樹冠に対して、1 つの樹木頂点が抽出されることが確認できた。
【０２６６】
隣接する樹木頂点を除去した結果を表５に示す。
【表５】
【０２６７】
H4-1、H4-2、H5-1 以外のプロットでは抽出精度に違いはないが、特にH4-1、H4-2 では樹頂点の抽出精度が向上した。現地で確認した樹冠形状とオルソ画像を勘案すると、H4 プロットは樹冠表層が他のプロットに比べて凹凸が明瞭ではなく、複雑な形状をしていると思われる。このため、1 つの樹冠に対して隣接する２つの樹冠上部が抽出され、その結果、１つの樹冠に対して隣接する２つの樹頂点が抽出されてしまったと考えられる。
【０２６８】
＜実施の形態３＞
前述の実施の形態１及び実施の形態２のＤＳＭは、 樹冠形状明確化を行うために樹冠内部の標高データを除去している。
【０２６９】
これは、樹冠を透過するパルスや航空機から斜めにレーザーが照射されることで樹木間の隙間から地盤に到達するパルスが存在するため、従来のDSM には樹冠内部の標高データ（樹冠表層の標高データではなく、樹冠等を透過して取得された標高データ）も含まれている（図４７参照）。このDSM を用いて樹高データ(DHM)を算出すると、樹冠形状が不明瞭になってしまう傾向にある。
【０２７０】
従って、樹木頂点抽出精度を向上させるためには、樹冠内部の標高データを除去し、樹冠形状を明確化したDSM を作成する必要がある。
【０２７１】
手動により樹冠内部の標高データを除去し、樹冠形状を明確化する手法もある。しかし、広範囲にわたる森林でのデータの処理には時間と労力がかかり、かつ、人為的な要因による再現性の問題が生じる可能性がある。
【０２７２】
先述したとおり、従来のDSM では地盤面の標高データも含まれているため、計測データは図４８（ａ）のようになっている。赤点で示す樹冠表層の標高データと青点で示す地盤面の標高データではその分布域に違いがあるので、この差を用いて地盤面に存在する標高データを除去することができる。
【０２７３】
反射強度は反射する物体により、その値が異なり、樹木などの植物の反射強度は大きいが、道路などの人工物や地盤などの土の反射強度は小さいという特徴を持つ。また、樹冠などを透過することで、反射強度は減衰する。この特徴を活かし、標高データだけでは除去できない樹冠内部に存在する標高データを除去する。
【０２７４】
図４８（ｂ）に示す青点は樹冠内部に存在する標高データであり、標高データだけでは除去できないデータである。赤点で示す標高データは樹冠表層に存在するため、反射強度は大きい。一方、青点で示す樹冠内部の標高データは樹冠を透過し、減衰しているために反射強度は小さくなる。この反射強度の違いを用いて樹冠内部のデータを除去することができる。
【０２７５】
以上より、自動処理による樹冠形状明確化に十分な成果が得られたので、本手法による樹冠形状明確化DSM を作成し、以後の解析（樹木頂点の抽出等）を行うのが好ましい。
【０２７６】
なお、本実施の形態１および実施の形態２は下記に説明するように利用されることもある。
【０２７７】
（一元管理化）
本実施の形態２により高精度に林分パラメータが取得できることがわかった。
【０２７８】
任意の林分に対して、「平均樹高」「立木密度」「推定樹冠面積」といった林分パラメータが得られる。さらにこのデータを用いることで間伐林の選定に役立つ「収量比数」や「適正立木密度」が得られる。また、レーザー計測データからは「林分傾斜角」や「レーザーパルスの地盤到達率」といった情報も得られる。データの一元管理化した例を以下に示す。
【０２７９】
[１］林分に対して複数の情報を付与できる。
【０２８０】
[２］上記の取得情報を取得情報の位置や地図等の情報と共に可視的に管理できる。
【０２８１】
[３]データ管理・更新が容易で、データの信頼性が高められる。
【０２８２】
[４]一元データ管理により、データの比較・解析・評価が容易である。
【０２８３】
[５]他の情報（森林簿）と併せてデータを集約し、その情報も一元管理できる。
【０２８４】
以下に収量比数、森林整備ランクについて述べる。
【０２８５】
1）収量比数
収量比数は密度管理図をもとに算出するもので、最多密度における幹材積に対する現存幹材積の割合である。具体的には樹高と立木密度から算出される。梨のスギ、ヒノキ林の密度管理図に描かれている曲線式を次に示す。
【０２８６】
スギ林
V = (0.07155968H-1.373859 + 5062H-2.869785 / N )-1 式（F）
logNRF = 5.370947 −1.495926logH 式（G）
ヒノキ林
V = (0.035147H-1.080773 + 4711.2H-2.922894 / N )-1 式（H）
logNRF = 5.7384 −1.8482121logH 式（I）
NRF：最多密度におけるha あたりの本数（本/ha）である。
【０２８７】
本研究ではプロット内に含まれる樹木本数（現地調査結果、樹冠未到達の中下層木も含む）と抽出本数（開発手法によるもので、表層木のみとなる）とをもとにした立木密度を用いて収量比数を算出した。その結果を図４９に示す。この図からプロット間の相対的な大小関係は大きく変わらないため、間伐の優先度を選定する情報として表層木のみの収量比数を用いることが可能とわかる。
【０２８８】
また、収量比数を0.7 として適正立木密度を算出した。この数値は「中庸仕立
て」と呼ばれる一般的な間伐の指標である。
【０２８９】
森林整備ランク
森林の生育状況だけでなく、樹種や森林の生育している状況にも注目して森林整備の優先度を評価する指標が森林整備ランクである。これは「森林管理ランク（収量比数）」「森林特性ランク（樹種や林齢）」「国土保全ランク（傾斜や保全対象）」を求め、総合的に森林整備の優先度を把握することができる。
【０２９０】
森林整備ランクを算出した結果を図５０，５１に示す。このことから、スギ林においては本指標により劣勢木の潜在状況も加味した森林整備計画が可能と考える。
【０２９１】
また、取得した情報を林分毎に整備し、GIS データとして一元管理する例を図５２ に示す。
【０２９２】
これにより、森林整備ランクをオルソ画像情報に載せて把握したり、プロットを選択することで図中右上のような林分の詳細データを把握したりできる。成可能であり、下記の様な利用が考えられる。
【０２９３】
[１] 森林の状態や森林所有者の把握が容易
[２] 優先度を考慮した効率的な森林整備計画
[３] 情報共有・公開社会に適した森林情報のデータ整備・運用
詳細な現地調査を実施することなく、簡易・広範囲
に森林整備の優先度が把握できることは、林業経営の観点から求められているコスト削減や要間伐林の選定が急務となっている現在の社会情勢に貢献するものと考える。また、林分パラメータだけでなく、レーザー計測データのオルソ画像や既存森林簿等も情報として加えることで、さらにGIS 一元管理の利便性が高まり、公開・共有情報化社会に貢献するものと考える。
【０２９４】
さらに、災害等により、地面が削りとられた場合に、本手法により予めこの個所の付近のそう本数を得ておいて、再びこの個所の付近の本数を求めて、その差を求めると、どの程度の木が川に流れたかが分る。
【０２９５】
（樹冠面積の推定）
樹冠面積推定は樹冠形状指数Ｅｅｉをもとに、watershed アルゴリズムで行っている。実施の形態１の樹冠面積推定手法では、あらかじめ樹冠サイズを設定し、樹冠面積を推定していたため、疎密度にばらつきのある林分においては十分な樹冠面積の推定ができていなかった。そこで、アルゴリズムの見直しを行い、推定精度の向上を試みた。
【０２９６】
樹冠面積推定手法の改良は実施の形態１と同様にwatershed アルゴリズムを用いて行った。実施の形態１では樹木頂点から周囲を検索して集水域を推定していたが、本利用例では閾値以上の任意の点がどの樹木頂点の集水域に属しているかという考え方に基づいて樹冠面積を推定している。
【０２９７】
本利用例では、樹冠形状指数Eeiに設定する閾値以上を対象に樹冠面積を推定するため、最適な閾値を設定する必要がある。そこで、最適な閾値を設定するために、樹高データと樹冠形状指数の断面形状をもとに最適値を決める。
【０２９８】
樹高（DHM）3m 以上を森林領域と定義したことから、これに基づき、樹高3m 以上の領域と樹冠形状指数のある値以上の領域とがほぼ合致するとき、この樹冠形状指数Ｅｅｉの値を閾値とする。図３３に樹高と樹冠形状指数の2 値化画像（背景にオルソ画像を表示）を示している。
【０２９９】
図３３で示したように樹高3m以上の領域に相当する樹冠形状指数は-60 であったため、樹冠形状指数Eeiの閾値を-60 として樹冠面積を推定する。
【０３００】
樹冠面積推定結果と精度改良した手法を用いて推定した樹冠面積を実施の形態１で作成した推定樹冠面積と比較した結果を図４２に示す。図を見ると実施の形態１の推定面積に比べて、樹冠の大きさに応じた推定面積となっている。また、目視判読の樹冠面積と同様の大きさであり、推定手法が改善されていることがわかる。
【０３０１】
本利用例をヒノキ林に適用した結果の画像を図４３に示し、検証結果を表６、表７、表８に示す。
【表６】
【表７】
【表８】
【０３０２】
表６から樹木の抽出数の精度は目視判読と比較してH4 林分で十分な精度ではなく（密林では目視判読結果は過小抽出となる）、また、被覆面積率（樹冠面積合計の抽出割合）は目視判読に比べて、全体的に過大であった。しかし、調査林分の状況やfirst パルスの地盤到達率からヒノキ林分は鬱閉していると考えられ、目視判読の樹冠面積が過小と考えられる（目視判読の樹冠面積は樹冠の大きさ、ヒノキの樹冠形状、判読写真の解像度の問題から、鬱閉状態を表すような判読面積に上方修正することは困難）。このため、林分面積をもとに樹冠面積の抽出精度を検証した。その結果、林分面積に対して推定樹冠面積の抽出率は概ね100%となり、推定樹冠が鬱閉状態を表していることがわかる。
【０３０３】
表７には単樹冠面積の比較を示したが、密林では目視判読が過小抽出となることや推定樹冠面積が目視判読に比べ過大であることから十分な精度が得られていない。
【０３０４】
また、目視判読の抽出精度を利用して林分内に含まれる実樹木本数を推定した結果及び、推定実樹木本数と開発手法による抽出樹木本数との比較結果を表８に示す。目視判読の抽出精度から推定する実樹木本数と開発手法による抽出樹木本数との割合は約9〜11 割となり、林分全域で高精度に樹木本数が抽出できていると考えられる。
【０３０５】
さらに、このことと推定樹冠が鬱閉状態を再現できていたことを考えると表７ の推定単樹冠面積の再現性は高い。
【図面の簡単な説明】
【０３０６】
【図１】本発明にいたる経緯を説明の説明図である。
【図２】本発明にいたる経緯を説明する説明図である。
【図３】本発明にいたる経緯を説明する説明図である。
【図４】本実施の形態による尾根谷画像と樹冠形状指数画像との比較説明図である。
【図５】本実施の形態の樹木頂点認識装置の概略構成図である。
【図６】本実施の形態を説明する概略フローチャートである。
【図７】検索範囲の設定を説明する説明図である。
【図８】本実施の形態の樹冠形状指数による画像を説明する説明図である。
【図９】樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【図１０】樹冠形状指数算出部３２のアスペクト比を求める処理を説明するフローチャートである。
【図１１】検索範囲における樹木高グリッドデータの抽出を説明する説明図である。
【図１２】更新するアスペクト比と角度の関係を説明する説明図である。
【図１３】アスペクト比の更新処理を説明するフローチャートである。
【図１４】抽出したＤｐｐｉと更新したアスペクト比等の情報を説明する説明図である。
【図１５】ファイル３７の樹冠形状指数情報Ｇｉの説明図である。
【図１６】樹冠形状指数の算出のフローチャートである。
【図１７】樹冠領域に対しての樹冠部との関係を説明する説明図である。
【図１８】局所領域の平均値（μ）と標準偏差（σ）との関係を説明する説明図である。
【図１９】具体的な樹冠部の抽出処理を説明するフローチャートである。
【図２０】頂点候補Ｄｇｉの説明図である。
【図２１】樹冠の上部のデータＤｈｉの説明図である。
【図２２】樹冠部の説明図である。
【図２３】本実施の形態による疎林の頂点抽出結果を説明する説明図である。
【図２４】本実施の形態の樹冠形状指数画像と尾根谷度画像との違いを説明するための説明図である。
【図２５】本実施の形態の樹冠形状指数画像と尾根谷度画像との違いを説明するための説明図である。
【図２６】本実施の形態による疎林の樹冠部抽出結果を説明する説明図である。
【図２７】本実施の形態の具体的なハードウエア構成図である。
【図２８】実施の形態２における樹冠形状指数の平滑処理（局所領域）による樹木頂点の減少を説明する説明図である。
【図２９】実施の形態２による森林領域抽出結果を説明する説明図である。
【図３０】本実施の形態２の樹木頂点認識装置の概略構成図である。
【図３１】本実施の形態２を説明する概略フローチャートである。
【図３２】本実施の形態２のガウシアンフィルタと局所平均の樹木頂点抽出結果の比較を説明する説明図である。
【図３３】本実施の形態２の利用例（樹冠面積の推定）の森林領域の最適閾値の設定根拠を説明する説明図である。
【図３４】本実施の形態２による樹冠部域の抽出結果を説明する説明図である。
【図３５】隣接頂点の不具合を説明する説明図である。
【図３６】本実施の形態２のアスペクト比更新処理を説明するフローチャートである。
【図３７】本実施の形態２の地上・地下開度の置き換え変更を説明する説明図である。
【図３８】本実施の形態２の抽出したＤｐｐｉと更新したアスペクト比等の情報を説明する説明図である。
【図３９】実施の形態２のファイル３７の樹冠形状指数情報Ｇｉｉの説明図である
【図４０】局所平均と局所加重平均を説明する説明図である。
【図４１】本実施の形態２のガウシアンフィルタの説明図である。
【図４２】本実施の形態２の利用例（樹冠面積の推定）の樹冠面積推定結果を説明する説明図である。
【図４３】本実施の形態２の利用例（樹冠面積の推定）の樹冠面積推定結果を説明する説明図である。
【図４４】本実施の形態２の隣接頂点樹木前後の比較を説明する説明図である。
【図４５】本実施の形態２の隣接樹木頂点除去処理を説明するフローチャートである。
【図４６】本実施の形態２の隣接樹木頂点の除去を説明する説明図である。
【図４７】実施の形態３のＤＳＭの標高データの削除を説明する説明図である。
【図４８】実施の形態３の樹冠内部データを説明する説明図である。
【図４９】本実施の形態２の利用例（収量比較）を説明する説明図である。
【図５０】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【図５１】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【図５２】本実施の形態２の利用例（森林整備ランク）を説明する説明図である。
【符号の説明】
【０３０７】
１０ データベース
１１ データベース
１３ 表示部
１４ 樹木高算出部
１７ 尾根谷度算出部
２１ 樹木頂点候補抽出部
２３ 樹冠部抽出部
２５ 第１の樹木頂点決定部
３０ 出力部
３２ 第１の樹冠形状指数算出部
３５ 閾値算出部
４１ 第２の樹冠形状指数算出部
４１ａ 第２のアスペクト比更新処理
４１ｂ 第２の樹冠形状計算処理
４１ｃ ガウシアンフィルタ
４１ｄ 森林領域抽出処理
４３ 第２の樹木頂点決定部
【特許請求の範囲】
【請求項１】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識方法であって、
コンピュータが、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶手段に記憶するステップと、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出ステップと、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出ステップと、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出ステップと、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識ステップと
を行うことを特徴とする樹木頂点認識方法。
【請求項２】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記検索方向ｋｉｐ内での各樹木高グリッドデータＤｐｐｉの樹木高と距離からアスペクト比を順次求め、これらのアスペクト比の中の最大値及び最小値を抽出し、該最大値を所定の第１の値に更新して前記地上開度φ1´を求めると共に前記最小値を所定の第２の値に更新して前記地下開度φ2´を求めて前記樹冠形状指数ＥＧｉを算出させることを特徴とする請求項１記載の樹木頂点認識方法。
【請求項３】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記所定エリア内の前記樹木高グリッドデータＤｄｉを順次設定するＡステップと、
該樹木高グリッドデータＤｄｉが設定される毎に、この樹木高グリッドデータに前記検索範囲Ｍｉを設定するＢステップと、
前記樹木高グリッドデータＤｄｉが設定される毎に、検索方向ｋｉｐを求めて設定するＣステップと、
前記検索方向ｋｉｐが設定される毎に、該設定された検索方向ｋｉｐに存在する前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉを順次設定し、水平軸上の距離に対しての垂直方向のアスペクト比を求め、求められたアスペクト比から最大および最小のアスペクト比を抽出するＤステップと、
前記Ｄステップで抽出された最大および最小のアスペクト比を基に、該最大および最小のアスペクト比を、前記第１の所定の値の針葉樹用の新たなアスペクト比及び前記第２の所定の値の針葉樹用の新たなアスペクト比に更新し、該更新した新たなアスペクト比で前記樹冠形状指数ＥＧｉを算出させるＥステップと
を行うことを特徴とする請求項２記載の樹木頂点認識方法。
【請求項４】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記樹木高グリッドデータＤｄｉに対して全ての検索方向ｋｉｐでの隣接するグリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスのアスペクト比になっているときは、前記Ｅステップで求められたアスペクト比を前記針葉樹の樹木頂点と認識しやすくなるように新たなアスペクト比にするＦステップと
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項５】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記アスペクト比の最大値及び最小値から前記地上開度φ１´及び地下開度φ２´を順次求めるＧステップと、
前記地上開度φ１´及び地下開度φ２´から樹冠形状指数φ３´を順次求めるＨステップと、
全ての前記検索方向ｋｉｐの前記樹冠形状指数φ３´が求められたときに、これらを平均化して、さらに局所加重平均化し、これを前記樹冠形状指数ＥＧｉとして記憶するＩステップと
を行うことを特徴とする請求項１記載の樹木頂点認識方法。
【請求項６】
前記コンピュータは、
前記樹冠部抽出ステップが、
前記樹冠形状指数ＥＧｉが記憶される毎に、該樹冠形状指数ＥＧｉを読み込むと共に、前記記憶手段の前記局所領域を読み込むＪステップと、
前記局所領域の平均値及び標準偏差の和を前記樹冠部の抽出のための閾値Ｈｇｉとして求めるＫステップと、
前記局所領域の中央の樹冠形状指数ＥＧｉが閾値以上の場合は、該中央の値を樹冠部とする値に置き換えて表示させるＬステップと
を行うことを特徴とする請求項１記載の樹木頂点認識方法。
【請求項７】
前記Ｄステップは、
前記樹木高グリッドデータＤｄｉ及び前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉが設定される毎に、垂直方向の差ΔＺｃｉ及び水平方向の距離ｄｉを求めるステップと、
前記差Ｚｃｉとその距離ｄｉとのアスペクト比を水平方向を基準にして順次求め、このアスペクト比を算出するステップと、
前記算出されたアスペクト比の中から最大と最小のアスペクト比を抽出するステップと、
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項８】
前記Ｅステップは、
前記Ｄステップで求められたアスペクト比の最大値又は最小値が、水平方向を基準にして前記垂直方向で、
前記アスペクト比の最大値又は最小値が−２．７４７４〜０（０を含まない）の範囲に相当する場合は、最大値又は最小値を、−２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が０〜＋２．７４７４に相当する場合は、最大値又は最小値を、＋２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が−２．７４７４以下に相当する場合は、最大値又は最小値を、＋２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が＋２．７４７４以上に相当する場合は、前記抽出された最大値又は最小値のままとするステップと
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項９】
前記Ｆステップは、
前記樹木高グリッドデータＤｐｐｉに対して全ての検索方向での隣接するグリッドデータとの比較で求められたアスペクト比が、前記水平軸を基準としてマイナスになっているときは、前記Ｅステップで求められたアスペクト比を１０００倍することを特徴とする請求項４記載の樹木頂点認識方法。
【請求項１０】
前記コンピュータは、
前記樹木頂点認識ステップが、
前記樹冠部毎に、該樹木部内に存在する最大の樹冠形状指数ＥＧｉの箇所を樹木頂点とするステップとを行うことを特徴とする請求項１乃至９のいずれかに記載の樹木頂点認識方法。
【請求項１１】
前記樹木頂点認識ステップは、
前記樹冠部毎に、ｎ×ｎのフィルタをかけて、該フィルタの中央に樹木頂点データが存在しかつ、該フィルタ内に他の樹木頂点データが存在したときに、両者の前記樹冠形状指数ＥＧｉを比較し、前記中央の樹冠形状指数ＥＧｉの方が小さい場合、前記中央の樹木頂点を消去することを特徴とする請求項１０記載の樹木頂点認識方法。
【請求項１２】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識装置であって、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶した記憶手段と、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で基づいて地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出手段と、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出手段と、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出手段と、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識手段と
を有することを特徴とする樹木頂点認識装置。
【請求項１３】
前記樹冠形状指数算出手段は、
前記検索方向ｋｉｐ内での各樹木高グリッドデータＤｐｐｉの樹木高と距離からアスペクト比を順次求め、これらのアスペクト比の中の最大値及び最小値を抽出し、該最大値を所定の第１の値に更新して前記地上開度φ１´を求めると共に前記最小値を所定の第２の値に更新して前記地下開度φ２´を求めることを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１４】
前記樹冠形状指数算出手段は、
前記所定エリア内の前記樹木高グリッドデータＤｄｉを順次設定するＡ手段と、
該樹木高グリッドデータＤｄｉが設定される毎に、この樹木高グリッドデータに前記検索範囲Ｍｉを設定するＢ手段と、
前記樹木高グリッドデータＤｄｉが設定される毎に、検索方向ｋｉｐを求めて設定するＣ手段と、
前記検索方向ｋｉｐが設定される毎に、該設定された検索方向ｋｉｐに存在する前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉを順次設定し、水平軸上の距離に対しての垂直方向のアスペクト比を求め、求められたアスペクト比から最大および最小のアスペクト比を抽出するＤ手段と、
前記Ｄ手段で抽出された最大および最小のアスペクト比を基に、該最大および最小のアスペクト比を、前記第１の所定の値の針葉樹用の新たなアスペクト比及び前記第２の所定の値の針葉樹用の新たなアスペクト比に更新し、該更新した新たなアスペクト比で前記樹冠形状指数を算出させるＥ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１５】
前記樹冠形状指数算出手段は、
前記樹木高グリッドデータＤｄｉに対して全ての検索方向ｋｉｐでの隣接するグリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスのアスペクト比になっているときは、前記Ｅステップで求められたアスペクト比を前記針葉樹の樹木頂点と認識しやすくなるように新たなアスペクト比にするＦ手段と
を有することを特徴とする請求項１３記載の樹木頂点認識装置。
【請求項１６】
前記樹冠形状指数算出手段は、
前記アスペクト比の最大値及び最小値から前記地上開度φ１´及び地下開度φ２´を順次求めるＧ手段と、
前記地上開度φ１´及び地下開度φ２´から樹冠形状指数φ３´を順次求めるＨ手段と、
全ての前記検索方向ｋｉｐの前記樹冠形状指数φ３´が求められたときに、これらを平均化して、さらに局所加重平均化し、これを前記樹冠形状指数ＥＧｉとして記憶するＩ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１７】
前記樹冠部抽出手段が、
前記樹冠形状指数ＥＧｉが記憶される毎に、該樹冠形状指数ＥＧｉを読み込むと共に、前記記憶手段の前記局所領域を読み込むＪ手段と、
前記局所領域の平均値及び標準偏差の和を前記樹冠部の抽出のための閾値Ｈｇｉとして求めるＫ手段と、
前記局所領域の中央の樹冠形状指数が閾値以上の場合は、該中央の値を樹冠部とする値に置き換えて表示させるＬ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１８】
前記Ｄ手段は、
前記樹木高グリッドデータＤｄｉ及び前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉが設定される毎に、垂直方向の差ΔＺｃｉ及び水平方向の距離ｄｉを求める手段と、
前記差Ｚｃｉとその距離ｄｉとのアスペクト比を水平方向を基準にして順次求め、このアスペクト比を算出する手段と、
前記算出されたアスペクト比の中から最大と最小のアスペクト比を抽出する手段と、
を有することを特徴とする請求項１７記載の樹木頂点認識装置。
【請求項１９】
前記Ｅ手段は、
前記Ｄ手段で求められたアスペクト比の最大値又は最小値が、水平方向を基準にして前記垂直方向で、
前記アスペクト比の最大値又は最小値が−２．７４７４〜０（０を含まない）の範囲に相当する場合は、最大値又は最小値を、−２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が０〜＋２．７４７４に相当する場合は、最大値又は最小値を、＋２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が−２．７４７４以下に相当する場合は、最大値又は最小値を、＋２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が＋２．７４７４以上に相当する場合は、前記抽出された最大値又は最小値のままとする手段と
を有することを特徴とする請求項１７記載の樹木頂点認識装置。
【請求項２０】
前記Ｆ手段は、
前記樹木高グリッドデータＤｐｐｉに対して全ての検索方向での隣接するグリッドデータとの比較で求められたアスペクト比が、前記水平軸を基準としてマイナスになっているときは、前記Ｅステップで求められたアスペクト比を１０００倍することを特徴とする請求項１９記載の樹木頂点認識装置。
【請求項２１】
前記樹木頂点認識手段は、
前記樹冠部毎に、該樹木部内に存在する最大の樹冠形状指数ＥＧｉの個所を樹木頂点とする手段
とを有することを請求項１２乃至２０のいずれかに記載の樹木頂点認識装置。
【請求項２２】
前記樹木頂点認識手段は、
前記樹冠部毎に、ｎ×ｎのフィルタをかけて、該フィルタの中央に樹木頂点データが存在しかつ、該フィルタ内に他の樹木頂点データが存在したときに、両者の前記樹冠形状指数ＥＧｉを比較し、前記中央の樹冠形状指数ＥＧｉの方が小さい場合、前記中央の樹木頂点を消去することを特徴とする請求項２１記載の樹木頂点認識装置。
【請求項２３】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識のプログラムであって、
コンピュータに、
樹冠形状指数を求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数から樹冠部を抽出するための局所領域を記憶手段に記憶する手段、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出手段、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出手段、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出手段、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識手段と
としての機能させるための樹木頂点認識のプログラム。
【請求項１】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識方法であって、
コンピュータが、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶手段に記憶するステップと、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出ステップと、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出ステップと、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出ステップと、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識ステップと
を行うことを特徴とする樹木頂点認識方法。
【請求項２】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記検索方向ｋｉｐ内での各樹木高グリッドデータＤｐｐｉの樹木高と距離からアスペクト比を順次求め、これらのアスペクト比の中の最大値及び最小値を抽出し、該最大値を所定の第１の値に更新して前記地上開度φ1´を求めると共に前記最小値を所定の第２の値に更新して前記地下開度φ2´を求めて前記樹冠形状指数ＥＧｉを算出させることを特徴とする請求項１記載の樹木頂点認識方法。
【請求項３】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記所定エリア内の前記樹木高グリッドデータＤｄｉを順次設定するＡステップと、
該樹木高グリッドデータＤｄｉが設定される毎に、この樹木高グリッドデータに前記検索範囲Ｍｉを設定するＢステップと、
前記樹木高グリッドデータＤｄｉが設定される毎に、検索方向ｋｉｐを求めて設定するＣステップと、
前記検索方向ｋｉｐが設定される毎に、該設定された検索方向ｋｉｐに存在する前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉを順次設定し、水平軸上の距離に対しての垂直方向のアスペクト比を求め、求められたアスペクト比から最大および最小のアスペクト比を抽出するＤステップと、
前記Ｄステップで抽出された最大および最小のアスペクト比を基に、該最大および最小のアスペクト比を、前記第１の所定の値の針葉樹用の新たなアスペクト比及び前記第２の所定の値の針葉樹用の新たなアスペクト比に更新し、該更新した新たなアスペクト比で前記樹冠形状指数ＥＧｉを算出させるＥステップと
を行うことを特徴とする請求項２記載の樹木頂点認識方法。
【請求項４】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記樹木高グリッドデータＤｄｉに対して全ての検索方向ｋｉｐでの隣接するグリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスのアスペクト比になっているときは、前記Ｅステップで求められたアスペクト比を前記針葉樹の樹木頂点と認識しやすくなるように新たなアスペクト比にするＦステップと
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項５】
前記コンピュータは、
前記樹冠形状指数算出ステップが、
前記アスペクト比の最大値及び最小値から前記地上開度φ１´及び地下開度φ２´を順次求めるＧステップと、
前記地上開度φ１´及び地下開度φ２´から樹冠形状指数φ３´を順次求めるＨステップと、
全ての前記検索方向ｋｉｐの前記樹冠形状指数φ３´が求められたときに、これらを平均化して、さらに局所加重平均化し、これを前記樹冠形状指数ＥＧｉとして記憶するＩステップと
を行うことを特徴とする請求項１記載の樹木頂点認識方法。
【請求項６】
前記コンピュータは、
前記樹冠部抽出ステップが、
前記樹冠形状指数ＥＧｉが記憶される毎に、該樹冠形状指数ＥＧｉを読み込むと共に、前記記憶手段の前記局所領域を読み込むＪステップと、
前記局所領域の平均値及び標準偏差の和を前記樹冠部の抽出のための閾値Ｈｇｉとして求めるＫステップと、
前記局所領域の中央の樹冠形状指数ＥＧｉが閾値以上の場合は、該中央の値を樹冠部とする値に置き換えて表示させるＬステップと
を行うことを特徴とする請求項１記載の樹木頂点認識方法。
【請求項７】
前記Ｄステップは、
前記樹木高グリッドデータＤｄｉ及び前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉが設定される毎に、垂直方向の差ΔＺｃｉ及び水平方向の距離ｄｉを求めるステップと、
前記差Ｚｃｉとその距離ｄｉとのアスペクト比を水平方向を基準にして順次求め、このアスペクト比を算出するステップと、
前記算出されたアスペクト比の中から最大と最小のアスペクト比を抽出するステップと、
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項８】
前記Ｅステップは、
前記Ｄステップで求められたアスペクト比の最大値又は最小値が、水平方向を基準にして前記垂直方向で、
前記アスペクト比の最大値又は最小値が−２．７４７４〜０（０を含まない）の範囲に相当する場合は、最大値又は最小値を、−２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が０〜＋２．７４７４に相当する場合は、最大値又は最小値を、＋２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が−２．７４７４以下に相当する場合は、最大値又は最小値を、＋２．７４７４に更新し、
前記アスペクト比の最大値又は最小値が＋２．７４７４以上に相当する場合は、前記抽出された最大値又は最小値のままとするステップと
を行うことを特徴とする請求項３記載の樹木頂点認識方法。
【請求項９】
前記Ｆステップは、
前記樹木高グリッドデータＤｐｐｉに対して全ての検索方向での隣接するグリッドデータとの比較で求められたアスペクト比が、前記水平軸を基準としてマイナスになっているときは、前記Ｅステップで求められたアスペクト比を１０００倍することを特徴とする請求項４記載の樹木頂点認識方法。
【請求項１０】
前記コンピュータは、
前記樹木頂点認識ステップが、
前記樹冠部毎に、該樹木部内に存在する最大の樹冠形状指数ＥＧｉの箇所を樹木頂点とするステップとを行うことを特徴とする請求項１乃至９のいずれかに記載の樹木頂点認識方法。
【請求項１１】
前記樹木頂点認識ステップは、
前記樹冠部毎に、ｎ×ｎのフィルタをかけて、該フィルタの中央に樹木頂点データが存在しかつ、該フィルタ内に他の樹木頂点データが存在したときに、両者の前記樹冠形状指数ＥＧｉを比較し、前記中央の樹冠形状指数ＥＧｉの方が小さい場合、前記中央の樹木頂点を消去することを特徴とする請求項１０記載の樹木頂点認識方法。
【請求項１２】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識装置であって、
樹冠形状指数ＥＧｉを求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数ＥＧｉから樹冠部を抽出するための局所領域を記憶した記憶手段と、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で基づいて地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出手段と、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出手段と、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出手段と、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識手段と
を有することを特徴とする樹木頂点認識装置。
【請求項１３】
前記樹冠形状指数算出手段は、
前記検索方向ｋｉｐ内での各樹木高グリッドデータＤｐｐｉの樹木高と距離からアスペクト比を順次求め、これらのアスペクト比の中の最大値及び最小値を抽出し、該最大値を所定の第１の値に更新して前記地上開度φ１´を求めると共に前記最小値を所定の第２の値に更新して前記地下開度φ２´を求めることを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１４】
前記樹冠形状指数算出手段は、
前記所定エリア内の前記樹木高グリッドデータＤｄｉを順次設定するＡ手段と、
該樹木高グリッドデータＤｄｉが設定される毎に、この樹木高グリッドデータに前記検索範囲Ｍｉを設定するＢ手段と、
前記樹木高グリッドデータＤｄｉが設定される毎に、検索方向ｋｉｐを求めて設定するＣ手段と、
前記検索方向ｋｉｐが設定される毎に、該設定された検索方向ｋｉｐに存在する前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉを順次設定し、水平軸上の距離に対しての垂直方向のアスペクト比を求め、求められたアスペクト比から最大および最小のアスペクト比を抽出するＤ手段と、
前記Ｄ手段で抽出された最大および最小のアスペクト比を基に、該最大および最小のアスペクト比を、前記第１の所定の値の針葉樹用の新たなアスペクト比及び前記第２の所定の値の針葉樹用の新たなアスペクト比に更新し、該更新した新たなアスペクト比で前記樹冠形状指数を算出させるＥ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１５】
前記樹冠形状指数算出手段は、
前記樹木高グリッドデータＤｄｉに対して全ての検索方向ｋｉｐでの隣接するグリッドデータのアスペクト比が、前記水平軸を基準として全てマイナスのアスペクト比になっているときは、前記Ｅステップで求められたアスペクト比を前記針葉樹の樹木頂点と認識しやすくなるように新たなアスペクト比にするＦ手段と
を有することを特徴とする請求項１３記載の樹木頂点認識装置。
【請求項１６】
前記樹冠形状指数算出手段は、
前記アスペクト比の最大値及び最小値から前記地上開度φ１´及び地下開度φ２´を順次求めるＧ手段と、
前記地上開度φ１´及び地下開度φ２´から樹冠形状指数φ３´を順次求めるＨ手段と、
全ての前記検索方向ｋｉｐの前記樹冠形状指数φ３´が求められたときに、これらを平均化して、さらに局所加重平均化し、これを前記樹冠形状指数ＥＧｉとして記憶するＩ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１７】
前記樹冠部抽出手段が、
前記樹冠形状指数ＥＧｉが記憶される毎に、該樹冠形状指数ＥＧｉを読み込むと共に、前記記憶手段の前記局所領域を読み込むＪ手段と、
前記局所領域の平均値及び標準偏差の和を前記樹冠部の抽出のための閾値Ｈｇｉとして求めるＫ手段と、
前記局所領域の中央の樹冠形状指数が閾値以上の場合は、該中央の値を樹冠部とする値に置き換えて表示させるＬ手段と
を有することを特徴とする請求項１２記載の樹木頂点認識装置。
【請求項１８】
前記Ｄ手段は、
前記樹木高グリッドデータＤｄｉ及び前記検索範囲Ｍｉ内の樹木高グリッドデータＤｐｐｉが設定される毎に、垂直方向の差ΔＺｃｉ及び水平方向の距離ｄｉを求める手段と、
前記差Ｚｃｉとその距離ｄｉとのアスペクト比を水平方向を基準にして順次求め、このアスペクト比を算出する手段と、
前記算出されたアスペクト比の中から最大と最小のアスペクト比を抽出する手段と、
を有することを特徴とする請求項１７記載の樹木頂点認識装置。
【請求項１９】
前記Ｅ手段は、
前記Ｄ手段で求められたアスペクト比の最大値又は最小値が、水平方向を基準にして前記垂直方向で、
前記アスペクト比の最大値又は最小値が−２．７４７４〜０（０を含まない）の範囲に相当する場合は、最大値又は最小値を、−２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が０〜＋２．７４７４に相当する場合は、最大値又は最小値を、＋２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が−２．７４７４以下に相当する場合は、最大値又は最小値を、＋２．７４７４に更新する手段と、
前記アスペクト比の最大値又は最小値が＋２．７４７４以上に相当する場合は、前記抽出された最大値又は最小値のままとする手段と
を有することを特徴とする請求項１７記載の樹木頂点認識装置。
【請求項２０】
前記Ｆ手段は、
前記樹木高グリッドデータＤｐｐｉに対して全ての検索方向での隣接するグリッドデータとの比較で求められたアスペクト比が、前記水平軸を基準としてマイナスになっているときは、前記Ｅステップで求められたアスペクト比を１０００倍することを特徴とする請求項１９記載の樹木頂点認識装置。
【請求項２１】
前記樹木頂点認識手段は、
前記樹冠部毎に、該樹木部内に存在する最大の樹冠形状指数ＥＧｉの個所を樹木頂点とする手段
とを有することを請求項１２乃至２０のいずれかに記載の樹木頂点認識装置。
【請求項２２】
前記樹木頂点認識手段は、
前記樹冠部毎に、ｎ×ｎのフィルタをかけて、該フィルタの中央に樹木頂点データが存在しかつ、該フィルタ内に他の樹木頂点データが存在したときに、両者の前記樹冠形状指数ＥＧｉを比較し、前記中央の樹冠形状指数ＥＧｉの方が小さい場合、前記中央の樹木頂点を消去することを特徴とする請求項２１記載の樹木頂点認識装置。
【請求項２３】
上空から地上を撮影した所定エリアのオルソフォト画像又は/及び所望のグリッドデータを画像化した画像を表示部の画面に表示する一方、前記所定エリア内の、地盤のＸ、Ｙ、Ｚ座標値及び樹木の表層のｘ、ｙ、ｚ座標値との差を求め、この差Ｚｃｉ及びＸｃｉ、Ｙｃｉ座標を樹木高グリッドデータＤｄｉとして得て、この樹木高グリッドデータＤｄｉに基づいて樹木の頂点を認識する樹木頂点認識のプログラムであって、
コンピュータに、
樹冠形状指数を求めるための前記グリッドデータの検索範囲Ｍｉ、前記樹木高グリッドデータＤｄｉ、前記樹冠形状指数から樹冠部を抽出するための局所領域を記憶手段に記憶する手段、
前記画像に対応する各樹木高グリッドデータＤｄｉを中心とした検索範囲Ｍｉ毎に、検索方向ｋｉｐでの各樹木高グリッドデータＤｐｐｉとの樹木高及び距離から求めたアスペクト比の最大値、最小値で地上開度φ１´、地下開度φ２´を求め、これらの開度に基づいて樹冠形状指数ＥＧｉを算出して前記記憶手段に記憶する樹冠形状指数算出手段、
前記樹木高グリッドデータＤｄｉに対して前記局所領域を定め、この局所領域毎に前記樹冠形状指数ＥＧｉの平均値と標準偏差とを求め、これらの和を閾値として求める閾値算出手段、
前記局所領域の前記閾値が求められる毎に、その樹冠形状指数ＥＧｉが前記閾値以上かどうかを判断し、閾値を超えているときに、樹冠部として抽出する樹冠部抽出手段、
前記樹冠部毎に、該樹木部内に存在する最適な樹木頂点を認識して、該認識した頂点を所定の形状にして前記表示部に表示又は前記頂点の座標値を出力する樹木頂点認識手段と
としての機能させるための樹木頂点認識のプログラム。
【図１】
【図２】
【図３】
【図５】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２７】
【図３０】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【図４１】
【図４５】
【図４６】
【図４８】
【図４】
【図６】
【図７】
【図８】
【図２３】
【図２４】
【図２５】
【図２６】
【図２８】
【図２９】
【図３１】
【図３２】
【図３３】
【図３４】
【図３５】
【図４２】
【図４３】
【図４４】
【図４７】
【図４９】
【図５０】
【図５１】
【図５２】
【図２】
【図３】
【図５】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２７】
【図３０】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【図４１】
【図４５】
【図４６】
【図４８】
【図４】
【図６】
【図７】
【図８】
【図２３】
【図２４】
【図２５】
【図２６】
【図２８】
【図２９】
【図３１】
【図３２】
【図３３】
【図３４】
【図３５】
【図４２】
【図４３】
【図４４】
【図４７】
【図４９】
【図５０】
【図５１】
【図５２】
【公開番号】特開２００９−２２２７８（Ｐ２００９−２２２７８Ａ）
【公開日】平成２１年２月５日（２００９．２．５）
【国際特許分類】
【出願番号】特願２００８−１６１０４７（Ｐ２００８−１６１０４７）
【出願日】平成２０年６月１９日（２００８．６．１９）
【出願人】（０００００３６８７）東京電力株式会社 (2,580)
【出願人】（５９１０７４１６１）アジア航測株式会社 (48)
【Ｆターム（参考）】
【公開日】平成２１年２月５日（２００９．２．５）
【国際特許分類】
【出願日】平成２０年６月１９日（２００８．６．１９）
【出願人】（０００００３６８７）東京電力株式会社 (2,580)
【出願人】（５９１０７４１６１）アジア航測株式会社 (48)
【Ｆターム（参考）】
[ Back to top ]