Silviase/jgov_v1.4_full_ann huggingface dataset:
# D-1006 熱帯林のREDDにおける生物多様性保護コベネフィットの最大化に関する研究

(1) リモートセンシングによる森林の3次元構造とその変化の把握手法の開発

（独）森林総合研究所

温暖化忒忑推進拋點溫暖化忒応推進室 平田 泰雅

森林管理研究领域 资源解析研究室 薔藤 英樹

<研究協力者>

森林管理研究领域 资源解析研究室 松浦 俊也

平成22～24年度累計予算額：18,322千円

（うち、平成24年度予算額：2,968千円）

予算額は、間接経費を含む。

## [要旨]

施業履歴に応じて森林劣化状況に違いが見られる熱帯生産林が広がるマレーシア・サバ州デラマコット、タンクラップおよびピナンガ森林管理区を対象に、衛星LiDARと高分解能衛星画像のオブジェクトベース画像解析および統計モデルを用いて、森林の3次元構造や地上部バイオマス量を推定するモデルを構築した。まず、森林の3次元構造を捉えるために、2003年から2009年に取得された衛星LiDARデータのうち、雲などの影響が少ないことを前処理により確認できた2008年の観測地点（直径70m）28ヶ所にGPSを用いて到達し、直径40mの円形プロットを設定し、胸高直径と樹高の毎木調査を行った。次に、衛星LiDARの観測波形から最大樹高と地面からの反射強度の相対積算値が10%、25%、50%、75%、90%の高さを求め、地上調査で得られた森林の3次元構造との関係を解析し、地上部バイオマスを推定するモデルを構築した。とくに、地面からの反射強度の相対積算値が50%の高さが地上部バイオマスや森林劣化状況と相関が強いことが明らかになった。さらに、森林劣化状況の違いを面的に捉えるために、高分解能衛星画像IKONOSを用いたオブジェクトベース画像解析により林冠構造の類似したエリアをゾーニングし、各バンドの平均値や標準偏差、林冠のテクスチャおよび地形条件等の様々な特徴量をエリアごとに算出し、一般化線形モデル等を用いて地上部バイオマスを推定するモデルを構築した。高分解能衛星画像は、一般にスペクトル分解能が低い（バンド数が少ない）弱点があるが、林冠のテクスチャや地形条件などのエリア毎の特徴量を加えることで、施業履歴の違いによる地上部バイオマスの相違を推定することが可能であった。さらに、このバイオマス推定結果を教師データとして2005年撮影のSPOT 5衛星画像のうちバイオマスが大きく変化した場所を抽出し、扱伐によって森林劣化が進んだ場所を検出できた。

## [キーワード]

衛星LiDAR、高分解能衛星、3次元構造、地上部バイオマス、熱帯低地林
### 1. はじめに

気候変動枠組条約における京都メカニズム後の新たな枠組みとして、REDD（Reducing Emissions from Deforestation and forest Degradation in developing countries: 発展途上国での森林減少・劣化に伴う温室効果ガス排出削減）が検討されている $ ^{1} $ 。REDD実施においては、国あるいは準国レベルで森林減少と森林劣化を評価する必要上、人工衛星データに基づく広域での森林炭素蓄積把握の技術開発が急務となっている $ ^{2} $ 。日本政府もその技術開発において国際ワークショップの開催や研究支援を通じた国際的貢献を行っているが、高バイオマスを有する熱帯林での森林劣化を高い精度で把握する技術は未確立である。

東南アジアの熱帯低地の生産林では、扱伐林施業により森林減少・劣化が進行している。扱伐林では、高木の伐採に加え、搬出のために多数の作業道が林内に作られ、これらの影響で森林の林冠構造が変化する³）。劣化した森林への再造林もなされつつあるが、衛星画像を用いて林冠構造の違いにもとづき林相を区分できれば、優先順位をつけた効率的な森林の再生を行う上で有効と考えられる。REDDと生物多様性のコベネフィットを両立させるためには、森林減少と森林劣化の両方の状況を捉える手法が求められる。ここで、樹冠よりも小さい1～4m程度の解像度をもつ高分解能衛星画像は林相把握に利用でき、オブジェクトベース画像分類による林冠構造把握にもとづく林相区分が各地の人工林や天然林でなされつつある⁴）。様々な施業履歴をもつ熱帯低地の生産林への高分解能衛星画像解析の適用可能性を調べることは、これらの地域における森林減少・劣化の状況の面的な把握において重要である。

また衛星LiDAR（Light Detection and Ranging）は、衛星から照射された近赤外のレーザー光が林冠表面から地表面までの各階層で反射して返ってくるまでの時間と反射強度を計測する技術であり、これにより森林の3次元構造を推定することが可能である $ ^{5} $ 。とりわけ、各階層の林冠からの反射強度を観測できることから（図(1)－1）、森林の高さの情報のみならず、森林劣化の度合いを評価することが期待される。

図(1)−1 衛星LiDARによる森林の3次元計測
## 2 ．研究開発目的

10 万欄造の把握手法を開発する。また、高分解能衛星から得られる林冠の面的構造と組み合わせることにより、リモートセンシングによる地上部バイオマス推定の高度化を図る。衛星 LiDAR は、衛星から照射された近赤外のレーザー光が林冠表面から地表面までの各階層で反射して返ってくるまでの時間と反射強度を計測する技術であり、これにより森林の3次元構造を推定することが可能である。また、現在高分解能衛星は地上分解能が1 m以下で林冠を構成する各樹冠を観測することができる。そこで、森林の3次元構造の異なるエリアに調査プロットを設置して3次元構造に関する調査を行い、衛星 LiDAR データにおける各階層からの反射強度との関係を明らかにする。さらに高分解衛星データから得られる樹冠の面積と対応する立木のバイオマスとの関係を相対成長を仮定することにより推定するとともに、林冠以外の階層のバイオマスを衛星 LiDAR から推定することにより、森林全体のバイオマス推定精度を向上させる手法を開発する。

まず、対象とする熱帯林管理区の林分特性を明らかにするため、高分解能衛星データを用いてオブジェクトベース画像分類により、林冠の構造の類似したエリアをゾーニングする。さらに、衛星 LiDAR データを整備すると共に、検証用のデータを得るため、衛星 LiDAR の観測点（フットプリント）の中心点の位置を GPS を用いて同定し、現地調査を行う。現地調査ではプロットを設定し、樹高や胸高直径、樹種などの森林の3次元構造を調査する。次に、衛星 LiDAR の各観測点でのレーザー光の各階層からの反射強度の波形を解析し、現地調査の結果との比較から、森林の3次元構造との関係を明らかにする。また、高分解能衛星データのゾーニング結果と地上での毎木調査の結果を統計解析し、林冠からみた林分構造とバイオマスの推定手法を開発する。最後に、高分解能衛星データから得られる林冠の面的構造と衛星 LiDAR データから得られる森林の3次元構造に関する情報から、リモートセンシングによるバイオマス推定の高度化を図る。

### 3. 研究開発方法

本研究のフローを図(1)-2に示す。

本研究では、衛星LiDARを用いた熱帯低地林の3次元構造の把握手法を開発し、衛星LiDARと高分解能衛星データを複合的に用いることで森林バイオマスの推定手法を高度化した。また、時系列で推定したバイオマスを比較することによりその変化をとらえた。衛星LiDARデータとしてはICESat（Ice, Cloud, and Land Elevation Satellite）衛星GLAS（Geoscience Laser Altimeter System）のデータを用いた。また、高分解能衛星としては、2009年のIKONOS衛星データおよび2003年のSPOT衛星5号のデータを用いた。

まず、衛星LiDARを用いて森林の3次元構造を把握するために、衛星で観測された時間に沿った反射強度を表す波形と現地調査から得られた林分の情報とを比較した。森林の3次元構造がこの波形に与える影響を明らかにし、波形の情報から林分因子を推定する手法を開発した。

次に、林分を類型化し、地上部バイオマス量の類似した林相区分から森林のバイオマスを推定するため、高分解能衛星データに対してオブジェクトベース画像分類を行った。オブジェクトベース画像分類においては、スペクトル反射特性に加え、林冠表面の凹凸が作り出すテクスチャや地形情報などを加えることにより、地上部バイオマス推定の高度化を図った。

最後に、高分解能衛星データ解析で得られた地上部バイオマス推定結果を教師データとして過
去の衛星データに対してバイオマスに関する画像分類を行い、過去の地上部バイオマス分布図を作成した。両者を比較することによりバイオマス変化を推定した。

図(1)-2 本研究のフロー

以下に、本研究で開発した手法について、詳細に述べる。

## （1） 研究对象地

熱帯低地の生産林が広がるマレーシア・サバ州中央部において、サバ州林業局が管理する異なる森林施業形態をもつ隣接する3ヶ所の森林管理区にまたがる地域（5°14-30°N、117°11-36°E）を対象とした（図(1)-3）。

1995年より低インパクト施業（Reduced Impact Logging）が進められているデラマコット森林管理区(C)では蓄積量の豊かな森林が広がっている（図(1)-4c） $ ^{6)} $ 。ここでは包括的伐採計画の策定が行われており、伐採後の河川への影響調査や、搬出重機による森林への影響を最小限に抑えるための搬出路の設定など、森林の健全性の維持に向けた配慮がなされている。

西隣するタンクラップ森林管理区(B)では、2000年頃まで従来型のインパクトの大きい択伐林施業が行われ、その後再植林が進むものの前者に比べて劣化した森林が広がっている（図(1)-4b） $ ^{7} $ 。

同区の南西に隣接するピナンガ森林管理区(A)は、2006～2007年に皆伐林施業が行われ（図(1)-4a）、その後再植林がなされている $ ^{7} $ 。これらの生産林地域の周囲には、1980年代頃からの大規模皆伐後に土地利用転換された広大なアブラヤシ園が広がっている。標高は数10mから300m程度で、東側にかけてやや標高が高くなる。
A

保護林

生產林

解析範囲(Ikonos画像)

B

注：(a)ピナンガ森林管理区、(b)タンクラップ森林管理区、(c)デラマコット森林管理区．各森林管理区の周囲には広大なアブラヤシプランテーションが広がる。

図(1)-3 対象地の位置図

(A) サバ州における生産林の分布。(B) 解析対象地周辺の Landsat TM 画像(2005 年)
(a) 皆伐地

(b) 抿伐林

(c) 低インパクト施業

図(1)-4 対象地の森林の様子

## （２）使用した衛星データ

本研究では、森林の3次元構造を捉えることが可能な衛星LiDARと単木樹冠を観測可能な高分解能衛星のデータを用いた。

まず、森林の3次元構造を波形解析から明らかにするための衛星LiDARデータとしてアメリカ航空宇宙局（NASA）が2003年1月12日に打ち上げたICESat衛星GLASデータを用いた。ICESat衛星GLASデータは、照射したレーザー光が地表面から反射してくるエネルギーを1ナノ秒（ns）単位で544回記録する波形型LiDARデータである。この衛星LiDARのフットプリントは直径が約70mであり、170m間隔でライン上に並んでいる。観測高は1nsから151nsまでは60cmに、152ns以降は15cmに相当する。データは、ICESatのホームページからリクエストし、アメリカ国立雪氷データセンター（NSIDC）によって処理されたものをダウンロードした。ダウンロードしたデータはASCIIデータに変換処理して統計ソフトを用いて解析した。解析には2008年10月8日に取得されたデータを用いた。GLASデータは14のデータセットから構成されている。このうち反射の波形の再現にはGLA01のデータを、また観測位置の地理座標の取得にはGLA14のデータを用いた（図(1)-5）。なお、ICESat衛星GLASは観測機器の故障により、2009年10月11日に観測を終了している。

次に、高分解能衛星データから森林を林冠の状態の類似した領域ごとにゾーニングし、林冠からみた林分構造と地上部バイオマスを推定する手法を開発するため、2009年9月16日に撮影された高分解能衛星IKONOSの画像を用いた。IKONOSには1m解像度のパンクロ画像と4m解像度の赤(R)、緑(G)、青(B)、近赤外(IR)のマルチバンドの観測波長域がある（表(1)-1）。解析に先立ち、画像の精密幾何補正を行い、さらに、画像内に散在する雲や影の場所にマスクをかけて解析から除外した。次に、各バンドに加え、植生の活力を表す正規化植生指数

$$ NDVI=(IR-R)/(IR+R) $$

を算出した。

また、バイオマスの面的変化を時系列で把握するため、2003年7月28日に撮影された高分解能衛星SPOT-5の画像を用いた。SPOT-5には2.5m解像度のパンクロ画像と10m解像度の赤(R)、緑(G)、
赤外(IR、SWIR)のマルチバンドの観測波長域がある（表(1)-1）。SPOT画像の幾何補正については時系列での変化を抽出するために用いることから、先のIKONOS衛星画像と一致するようにIKONOS画像を参照画像として幾何補正を行った。

=GLA01
version=v5_6
type=LONG
i_rec_ndx	=	1067299551
i_UTCTime	=	308111110	475249
i_gla01_rectype	=	2
i_Spare1	=	0
i_fillnum	=	0	0	0	0	0	0	0
i_shot_ctr	=	9	10	11	12	13	14	15
16
i_statflags	=	66326527	66326527	66326527	66326527
66326527	66326527	66326527
i_gainSet1064	=	250	250	250	250
250
i_4nsPeakVal	=	32	32	32	32
32
i_8nsPeakVal	=	30	30	30	30
30
i_4nsBgMean	=	2936	2952	2947	2936
2936
i_4nsBgSDev	=	139	145	143	140
140
i_samp_pad	=	44	44	44	44
44
i_comp_type	=	1	1	1	1	1	1	1
i_rng_wf	=	27	29	29	29	30	30	30	32	31	28	29
28	28	29	30	28	29	30	32	32	31	28	29	32
29	28	29	29	28	28	28	31	28	28	29	28	29
31	29	29	31	29	27	32	30	26	29	27	28	30
28	28	29	30	30	30	28	31	32	29	28	28	30
31	30	28	29	28	28	29	29	28	30	31	30	28
29	31	30	30	32	31	28	31	30	28	29	31	29
27	28	30	30	31	31	29	31	28	29	29	28	30
30	30	30	28	29	30	30	30	31	29	29	28	29
29	28	27	29	28	29	29	28	29	30	28	29	30
28	29	31	30	27	29	29	28	29	29	28	29	30
28	29	29	29	32	29	28	28	31	30	27	28	29
29	29	29	29	29	29	28	31	30	27	27	28	29
28	29	29	29	29	28	28	31	30	27	27	28	29
28	29	29	29	29	28	28	31	30	27	27	28	29
28	29	29	29	29	28	28	31	30	27	27	28	29
<ecel>

図(1)-5 ダウンロードしたICESat衛星GLASGLA・GLA01データのASCIIへの変換

表(1)-1 IKONOS衛星とSPOT衛星の観測波長帯

バンド	IKONOS衛星	SPOT衛星5号
青	0.445-0.516 $ \mu $ m	-
緑	0.506-0.595 $ \mu $ m	0.50～0.59 $ \mu $ m
赤	0.632-0.698 $ \mu $ m	0.61～0.68 $ \mu $ m
近赤外	0.757-0.853 $ \mu $ m	0.78～0.89 $ \mu $ m
短波長赤外	-	1.58～1.75 $ \mu $ m
パンクロ	0.45-0.90 $ \mu $ m	0.48～0.71 $ \mu $ m
## （３）プロット設定と調査

衛星 LiDAR による観測は、衛星の軌道に沿ったライン上の観測となっている。軌道に沿って観測点（フットプリント）が 170m 間隔で並ぶ一方、軌道間の距離は数 10 km に及ぶため、本対象地ではタンクラップ森林管理区にのみ観測点（フットプリント）が存在した。そこで、衛星 LiDAR による森林の 3 次元構造解析用の参照データと、高分解能衛星データによる画像分類用のグランドトゥルースデータを取得するため、対象地域内にそれぞれの目的に応じた暫定プロットを設定した。衛星 LiDAR データの解析用に設定したプロットデータのうち IKONOS 画像の範囲に含まれるものについては画像分類にも用い、それぞれのプロットにおいて共通の項目について調査を行った。また、サバ州林業局森林研究所の厚意により、サバ州林業局森林研究所が対象地域で実施したプロット調査のデータも使用することができた。

まず、衛星 LiDAR データ解析用の暫定プロットを衛星 LiDAR のフットプリントと中心が一致するように設置した（図(1)-6）。各フットプリントの中心点には、予め LiDAR データから算出した座標に GPS のナビゲーション機能を使って到達し、現地でフットプリントの中心点を中心とした直径 40m の円形プロットを設定した。円形プロットの設定には、距離計測の可能な Vertex IV（Haglof 社製、スウェーデン）を用い、なるべく下層植生の刈り払いを少なくしたうえで、フットプリントの中心からプロット周囲を確定した。距離は水平距離とした。また、プロットの位置の正確性を高めるために、プロットの中心点において GPS を用いてプロット設定および毎木調査実施中に 1 時間以上の測位データの平均化を行った。この座標を GIS にインポートし、すべてのプロットが衛星 LiDAR のフットプリント内に設置されていることを確認した。胸高直径 10cm 以上のすべての樹木について、直径テープを用いて胸高直径（DBH）を、Vertex IV を用いて樹高（H）を測定した（図(1)-7）。樹種については、先駆樹種で非常に比重の軽い Macaranga かどうかについて野帳に記述した。

ここで得られたデータはプロットごとに表計算ソフトで整理・集計した（表(1)-2）。現地調査は、マレーシア・サバ大学国際森林学科の協力を得て、2010年11月、2011年2月、10月、2012年2月、2012年12月に実施した（図(1)-7）。最終的に、28か所の暫定プロットを衛星LiDARのフットプリント内に設置した。

高分解能衛星データの分類のためのグランドトゥルースデータとしては、2010～2012年に設定した円形プロット（半径20mで35ヶ所、半径15mで29ヶ所）と、サバ州林業局森林研究所およびデラマコット森林管理区が2007～2008年に設置した半径20m円形プロット（計69ヶ所）における毎木調査結果（DBHが10cm以上）を用い、このうち、IKONOS画像の範囲内に含まれ、かつ雲や影等の影響が少ないプロットを抽出した（図(1)-8）。衛星LiDARデータ解析用のプロットが全て扱伐林（タンクラップ森林管理区）に設定されているため、低インパクト施業が行われているデラマコット森林管理区（高バイオマス）や皆伐跡地のピナンガ森林管理区（低バイオマス）にもプロットを配置した。なお、プロット内に池が含まれた3プロットは解析から除外し、サバ州林業局森林研究所のプロットについては、本研究の共同研究者が現地確認できたプロットのみを用いた。また、同一の画像セグメントに含まれたプロット（2ヶ所）では、地上部バイオマス量の平均値を求めた。用いたプロット数は計41ヶ所である。DBHからの地上部バイオマスの算出には、熱帯雨林地域において作成されてきた複数の一般アロメトリ式を用いた（表(1)-3）。なお、本研究で新たに設置したプロットでは樹種を記録していないため、Chave et al式では材密度を一定
図(1)-6 衛星 LiDAR 解析用の暫定プロットの配置 衛星 LiDAR の観測点上に調査プロットを設定した。

図(1)-7 現地調査の様子（左：樹高の計測、右：胸高直径の測定）
表(1)-2 現地調査にて取得したプロットの概要

Plot No.	立木密度	地上部バイオマス(t/ha)	最小值	中位数	最大值	平均值	標準偏差	20m以上の本数
18-02	525	331	7.9	14.8	38.9	16.7	6.9	17
18-07	263	132	4.0	11.7	34.0	14.2	7.9	7
18-08	350	154	8.1	16.0	36.8	17.2	6.0	13
18-09	350	202	6.7	14.5	36.4	16.5	7.3	7
18-10	223	197	5.2	15.0	37.9	17.3	9.4	9
18-11	501	426	6.5	18.9	47.5	21.1	9.0	29
18-12	509	278	6.1	16.6	48.0	18.5	7.2	20
23-28	406	325	5.8	13.0	45.6	15.3	8.0	11
23-29	398	412	8.4	14.6	69.5	19.2	11.7	13
23-30	517	795	7.7	15.5	45.1	19.1	9.1	24
23-31	271	98	5.6	12.3	35.4	14.5	7.1	7
23-32	414	195	7.2	15.8	33.4	16.6	6.3	13
28-12	541	428	2.9	16.5	39.0	19.3	7.8	24
28-13	621	309	7.4	14.5	45.2	16.5	6.9	18
28-14	859	289	9.6	17.3	42.4	18.0	5.4	28
28-15	525	266	6.0	15.9	40.7	17.0	6.8	18
28-16	589	179	6.6	16.1	32.0	15.7	4.3	9
28-17	382	172	8.5	14.7	28.7	15.9	5.2	11
28-18	422	243	6.0	16.2	44.0	17.8	7.9	17
28-19	477	239	7.3	12.4	44.8	14.2	6.6	7
28-21	462	158	9.1	17.5	40.9	18.9	7.4	24
28-22	676	196	7.3	16.2	43.0	17.9	6.4	20
28-23	660	181	8.5	15.8	47.0	17.5	7.1	15
28-24	692	252	9.2	18.4	46.7	19.9	6.4	27
注：(A) ピナンガ森林管理区（皆伐地あり）、(B) タンクラップ森林管理区（2000年まで従来型の高インパクト扱伐後に再植林）、(C) デラマコット森林管理区（1995年から低インパクト施業(RIL)）。IKONOS ©

図(1)-8 高分解能衛星データ解析用の暫定プロットの配置

表(1)-3 一般アロメトリ式

Forest type	Equation	Reference
1 Asian moist forest (a)	$ \exp(-2.134 + 2.530 \times \ln(D)) $	Brown (1997)
2 Asian moist forest (b)	$ 42.69 - 12.8 \times D + 1.242 \times D^{2} $	Brown (1997)
3 Moist forest, $ p = 0.5 $	$ p \times \exp(-1.499 + 2.148 \times \ln(D) + 0.207 \times (\ln(D))^{2} - 0.02081 \times (\ln(D))^{3}) $	Chave et al. (2005)
4 Moist forest, Borneo	$ \exp(2.62 \times D - 2.3) $	Yamakura et al. (1986)
5 Degraded moist forest, Borneo	$ 0.0829 \times D^{2} + 2.43 $	Kenzo et al. (2009)

D: dbh, p: wood density

## （４）地上調査による樹高と地上バイオマスとの関係

衛星 LiDAR から得られる樹高に関する情報からバイオマスが推定可能であるかを検証するため、まず胸高直径を用いて地上部バイオマス量（AGB）を Brown の式（表(1)-3-1）より求めた。次にプロット内の立木密度（本/ha）を計算し、さらに樹高データからプロット内の最小値、第2四分位数、中位数、第3四分位数、最大値、平均値、標準偏差および20m以上の本数などの樹高関連変数を計算した（表(1)-2）。これらのデータから地上部バイオマス量を目的変数、樹高関連変数を説明変数とした重回帰分析を行い、樹高関連変数からの地上部バイオマス量推定可能性を調べた。また、樹高変数を用いてクラスタ分析を行い、得られたクラスタのバイオマスと樹高変数の関連を調べ、樹高変数を用いてプロットをバイオマスに応じたクラスに分類可能かどうかを調べた。

現地調査では、多数の樹木について胸高直径と樹高の両方を測定したが、一般に熱帯林で
の調査では胸高直径の測定しか行われていない。このため胸高直径（DBH）から樹高（H）を予測するモデルの検討を行った。今回は既存の DBH-H 関係を示すモデルとして、指数モデル

$$ \mathrm{H}=\mathbf{a}+\mathbf{b}*\log(\mathrm{DBH}) $$

と、非線形モデル⁸)

$$ \mathrm{H}=1/(1/(\mathrm{a^{*}}\mathrm{DBH})+1/\mathrm{b}) $$

を選定し、これらへ現地調査で測定したデータをあてはめモデルからパラメータを算出した。図(1)-9 に示すように、すべての樹木の DBH -H パターンと、低バイオマス林(-200Mg/ha)、中バイオマス林(200-300Mg/ha)、高バイオマス林(300Mg/ha-)では DBH -H パターンは異なっている。このため、あてはめに際しては、すべての樹木のデータを用いたあてはめだけでなく、森林劣化の程度に応じてプロットを低バイオマス、中バイオマス、高バイオマスに層化し、それぞれについてあてはめを行った。あてはめた結果については赤池の情報量基準（AIC）を用いて評価し、直径のみを測定したプロットから森林の3次元構造が再現可能であるかどうか調べた。

(a)

(b)

(c)

(d)

図(1)-9 直径と樹高の関係

(a) 低バイオマス林、(b) 中バイオマス林、(c) 高バイオマス林、(d) すべての樹木
## （５）衛星LiDARを用いた森林の３次元構造の計測

## 1 ）森林の林冠の状態と衛星LiDARデータとの関係

まず、ICESat衛星GLASデータから森林の3次元構造を把握するため、レーザー光の反射強度の波形データをGLA01データから抽出するプログラムを作成し、現地調査を実施した観測点（フットプリント）における反射強度の波形データを切り出した。さらに波形を可視化し、熱帯林の3次元構造をとらえる際のGLASデータの特徴を明らかにした。反射強度の情報にはノイズが含まれるが、この段階ではノイズ量の評価は行わず、可視化された波形のパターンだけを用いた。

地上では把握困難な林冠の状態について、高分解能衛星画像とフットプリントを重ね合わせ、高分解能衛星画像から得られた樹冠の状態と比較した。衛星LiDARの観測範囲であるタンクラップ森林管理区では択伐が行われているが、択伐の強さは必ずしも一様ではなく、非常に劣化した森林から比較的林冠が閉鎖した状態の森林まで、劣化の段階の異なる森林が見られる。そこで、劣化の状態に着目して、衛星LiDARの観測波形が劣化の状態によりどのように変化するのかを調べた。

## 2 ）最大樹高の推定

衛星LiDARによる観測において、衛星から照射されたレーザーは、フットプリント内の最大高の樹木の樹冠で反射を開始し、地面で最後に反射する $ ^{9)} $ （図(1)-10）。このため、衛星LiDARの反射の最初の信号と最後の信号の時間差から、林分の最大樹高が推定できると期待される。林分の最大樹高は、森林の遷移段階、立地環境、地上部バイオマスなどを推定するための有効なパラメータである。そこで、現地調査を実施した暫定プロットを含むフットプリントの衛星LiDARの観測波形データから、反射の信号の開始と終了を調べて最大樹高を推定した。

衛星LiDARのレーザーの反射強度の値にはノイズ情報が含まれているため、波形データからノイズ部分を推定し、ノイズの値を引いた値を各時間に観測された反射強度から減じた。この反射強度を観測時間に対してプロットし、反射の最後のピークの地点を地面とし、ノイズを引いた値が0またはマイナスでなくなる地点を樹木の最高地点とした（図(1)-10）。ここで推定した最大樹高を現地調査の暫定プロットで得られた最大樹高と比較した。

## 3 ）地上部バイオマスの推定

衛星LiDARで観測された反射強度の波形は、森林の林冠表面の構造によって形成される。すなわち、上空から観測可能な林冠の高さと階層構造により波形が決定される。森林の地上部バイオマスは、林冠を構成する樹木の高さや林冠の開空率、森林の階層構造と強い相関がある。衛星LiDARが観測する反射強度は、地面からのそれぞれの高さにおける上空から観測可能な林冠の面積と密接に関係していると考えられる。そこで、この反射強度の積算値を用いて地上部バイオマス推定モデルを作成した。

衛星LiDARの各観測点でのノイズを取り除いたエネルギーの反射強度について、地面からの最大樹高までの相対積算値が10%、25%、50%、75%、90%の高さ $ (h_{10}, h_{25}, h_{50}, h_{75}, h_{90}) $ を算出し（図(1)-10）、これらの変数を用いて地上部バイオマス（AGB）を推定するための次式の重回帰モデルを作成した。
$$ \mathrm{A G B}=\beta_{0}{h_{10}}^{\beta_{1}}{h_{25}}^{\beta_{2}}{h_{50}}^{\beta_{3}}{h_{75}}^{\beta_{4}}{h_{90}}^{\beta_{5}} $$

(1)-4 式は両辺の対数をとり、

$$ \ln\left(AGB\right)=ln\beta_{0}+\beta_{1}ln\ h_{10}+\beta_{2}ln\ h_{25}+\beta_{3}ln\ h_{50}+\beta_{4}ln\ h_{75}+\beta_{5}ln\ h_{90} $$

となる。このモデルから、どの高さでの反射が地上部バイオマスの推定に影響しているかを調べた。

図(1)-10 GLAS データの反射強度と地面からの最大樹高までの反射強度の相対積算值

## （６）高分解能衛星データを用いたバイオマス推定手法の開発

## 1 ）解析の流れ

解析方法の概要を図(1)-11 に示す。まず、対象地内に多数の毎木調査プロットを設定し、胸高直径(DBH)10cm 以上の全樹木について DBH を計測し、熱帯雨林地域において構築された、DBH のみにもとづく複数の一般アロメトリ式を用いてプロット毎の地上部バイオマスを算出した。次に、高分解能衛星画像 IKONOS と 30m DEM（SRTM:シャトルレーダートポグラフィックミッションにより作成されたもの）を用い、衛星画像を林冠状態が類似する場所に領域分割した上で、分割された領域ごとに、衛星画像の各バンドの特徴や地形条件を示す多数の特徴量を算出した。さらに、プロット毎の地上部バイオマスと、プロットが含まれる画像領域における各種特徴量との関係を、分類回帰木（CART: Classification and regression tree）や一般化線形モデル等の統計モデルを用いて解析し、地上部バイオマスの大小に関わりのある特徴量を調べた。そして、得られた回帰式を、調査プロットがない多数の画像領域に外挿し、対象地全体における地上部バイオマスの分布特徴を捉えた。以下に解析内容を詳述する。
図(1)-11 解析の流れ

## 2 ）IKONOS 衛星画像の前処理とオブジェクトベース画像解析

2009 年 9 月 16 日に撮影された高分解能衛星画像 IKONOS を用いた。IKONOS には 1m 解像度のパンクロ画像と 4m 解像度の赤(R)、緑(G)、青(B)、近赤外(IR)のマルチバンドの観測波長域がある。まず、衛星画像の RPC ファイルと SRTM 30m (30m DEM) を用いて衛星画像を幾何補正した。次に、マルチバンド（4 バンド）の画素値を用い、オブジェクトベース画像解析ソフトウェア cCognition Developer 8.7(Trimble 社)を用いて画像を領域分割した。領域分割には多重解像度領域分割法を用いた。この手法は、region growing アルゴリズム（領域拡張法）にもとづき画像内で近接する似た画素を同じ領域（セグメント）に分けていくもので、人間の判読結果に近い領域分割結果を得られることからオブジェクトベース画像解析において最もよく用いられる。eCognition では、いくつかのパラメータを設定することで解析目的に適した大きさや形状の領域分割が行える。本研究では、森林に関わる既往研究を参考に、カラー、シェープ、コンパクト性の各パラメータを 0.5 に設定した。ここで、カラーとシェープのパラメータは画素値と分割領域の形状の複雑さとの相対的な重要度を決めるものであり、両パラメータを合計すると 1 になる。森林においては形状よりも個々の画素値が重要となるが、高分解能衛星画像において画素値を重視しすぎると、樹冠の影やギャップに沿った複雑な形状のセグメントが生成されやすくなり、後述する毎木調査プロットが単一領域に収まらない場合が多くなる。そこで、カラーとシェープを同等に重み付けることで、個々の毎木調査プロットができるだけ一つのセグメントに含まれやすいようにした。また、生成される分割領域の面積に関わるスケールパラメータ(SP)を段階的に変化させて、領域分割結果を黙示判読で比較し、林冠構造の違いにより林相を良好に区分できた SP=50 にて領
域分割した。このときセグメントの平均面積は1.19 haとなり、後述する円形調査プロット(半径15mで0.07ha、半径20mで0.13ha)よりも十分大きいセグメントが生成された。SPを変えることで領域分割結果が変化する様子を図(1)-12に示す。なお、用いる画像の解像度に応じて適したSPの大きさは変化する。IKONOSのパンクロマティックバンド(解像度1m)を含めて解析した場合、同程度の面積の領域分割には、マルチバンド（解像度4m）を用いた場合に比べて4倍程度大きいSP=200を用いると、同程度の大きさ（平均面積0.83ha）の領域分割が行えた。なお、CARTによる解析ではパンクロマティックバンドを含めて解析したためSP=200の領域分割結果を用いたが、同バンドを含めることが解析結果の向上に寄与しないことが分かったため、一般化線形モデルによる解析ではマルチスペクトルバンドのみを用いてSP=50による領域分割結果を利用した。

IKONOS 画像内には多数の鱗状の雲とその影が散在しており、これらの場所は解析に利用できないため、マスクを作成して解析対象から除外した。ここで、高分解能衛星画像においては、個々の樹冠の影や林冠のギャップと雲の影とをピクセル単位で区別することが難しい。そこで、領域分割の結果を利用し、各領域内における青バンドと近赤外バンドの平均値について閾値を設定し、雲やヘイズ、裸地や道路などの明るい場所を青バンドにて、影や水域を近赤外バンドにて閾値を設定してマスクをかけた。さらに、マスク結果と元画像を重ねあわせて目視判読し、マスクの過不足がある場合はセグメントごとに追加または削除することで調整した。このようにして、個々の樹冠の影やギャップなどの小面積の特徴を残しつつ、雲やヘイズ、雲の影、裸地、水域などの十分な面積をもつセグメントのみにマスクをかけた。

A

B

図(1)-12 スケールパラメータ(SP)を変化させた場合の領域分割結果. (A) SP=30, (B) SP=50

注1：IKONOS©マルチスペクトル画像（空間解像度4m）を対象地の一部を切り出したもの．オレンジの線はオブジェクトベース画像解析による領域分割の境界線を示す。

注2：(a) ピナンガ森林管理区の2006年以降の皆伐・植林地、(b) タンクラップ森林管理区の従来型扱伐跡地、(c) タンクラップ森林管理区内の保存林
表(1)-4 說明変数

Variables	Equation
Mean and standard deviation (SD) of digital number (DN)
Blue, Green, Red, NIR
Mean and SD of vegetation index
NDVI	(NIR-Red)/(NIR+Red)
SVI	NIR/Red
GNDVI	(NIR-Green)/(NIR+Green)
TNDVI	(NDVI+0.5)0.5
Terrain feature derived from SRTM 30m
Mean_slope_position	Mean {elevation (center) - Mean elevation (7*7)}
Mean_gradient
Mean TWI (topographic wetness index)	Mean {ln (A/tan B)}
Area ratio of pixel-based ISODATA class within a segment 1)
shadow ratio
bareland_ratio
darker_vege_ratio
lighter vege ratio
Grey level co-occurrence matrix (GLCM) based texture variable
Angular second moment	∑i ∑j {p(i,j)}2
Contrast	∑n=0n2{∑i=1n∑j=1p(i,j)}
Correlation	∑i ∑j (ij)p(i,j)−μxμy
Dissimilarity	∑n=0n{∑i=1n∑j=1p(i,j)}
Entropy	− ∑i ∑j p(i,j)log(p(i,j))
Homogeneity	∑i ∑j 11+(i−j)2p(i,j)

15 classes --> 4 classes (shadow, bareland, darker vegetation, lighter vegetation)

## 3 ）セグメント毎の各種特徴量の算出

オブジェクトベース画像解析により生成したセグメント毎に算出した各種特徴量を表(1)-4に示す。まず、IKONOSマルチスペクトルの4つのバンド（青、緑、赤、近赤外）について、セグメント毎に平均と標準偏差（SD）を求めた。また、正規化植生指数（NDVI）等の各種植生指数についても、セグメント毎に平均と標準偏差を求めた。樹冠の影や林冠のギャップの割合と地上部バイオマスとの関わりを見るために、対象画像全体についてピクセルベースでのISODATA法
に基づく教師なし分類を行い、影や裸地の面積割合をセグメント毎に算出した。また、地上部バイオマスとテクスチャとの関係をみるために、グレーレベル同時生起行列（GLCM）を用いた各種テクスチャ情報を各バンドについてセグメント毎に算出した。その他、地形条件と地上部バイオマスとの関係をみるために、傾斜度、各ピクセルと周囲の7×7ピクセルの平均標高との差から求めた斜面相対位置、地形湿性指数（TWI）の3指標を30m DEMから算出し、セグメント毎の平均値を求めた。

次に、個々の毎木調査プロットと、オブジェクトベース画像解析で生成されたセグメントを重ねあわせ、単一のセグメントにプロットが完全に含まれる場合は、同プロットがそのセグメントを代表しているとみなして各種特徴量を対応づけた。部分的に複数のセグメントにまたがるプロットについては、円形プロット領域内に占める各セグメントの面積割合を求め、面積按分によって、各プロットとセグメントの特徴量を対応づけた。

## 4 ）分類回帰木（CART）による解析

ノンパラメトリックな解析手法のひとつである CART を用いて、IKONOS 画像と地上部バイオマスとの関係を解析した。CART では、各説明変数の二分割を繰り返し、目的変数ごとのばらつきが最小となるように分類基準が選択される。この手法は、欠損値をもつデータや非線形のデータも扱うことができ、得られた結果を解釈しやすいなどの利点をもつ。分類には統計ソフトウェア R.2.12.1 の mvpart パッケージを用いた。目的変数には、前述の直径 30m または 40m の円形の毎木調査プロットデータから Brown (a) によるアロメトリ式（表(1)-3 参照）により求めた地上部バイオマスの値をあらかじめ大中小に 3 区分したものを用いた（小：<150、中：150-300、大：300<、単位：t/ha）。一方、説明変数には、各円形プロット内における IKONOS 画像の各バンドの DN と NDVI それぞれの平均値と標準偏差を用いた。さらにクロスバリデーションにより、CART により生成された樹形図を樹形の複雑さと予測誤差の関係に基づき剪定した。

## 5 ）一般化線形モデルによる統計解析

地上部バイオマスと各種特徴量との対応関係の把握には一般化線形モデル（GLM）を用いた。GLMは線形回帰モデルを拡張したもので、リンク関数を変えることで様々な確率分布を仮定したモデルを扱うことができる。本項では、リンク関数にガウジアンのidentity（通常の線形回帰）、inverse（反比例）、log（ポアソン回帰）の3つを用いた。目的変数には各種アロメトリ式から算出した地上部バイオマスの整数値を用い、説明変数にはプロット毎の各種特徴量を用いた。ここで、ポアソン回帰は非負の整数値を目的変数とするモデルであるが、プロット間の地上部バイオマスの小数点以下での差異は無視できる程十分に小さいため、地上部バイオマスを整数值とすることで利用しうると考えた。

統計解析には R version 2.14.0 を用い、単回帰で有意(p<0.01)となる変数について重回帰にあてはめ、ステップワイズ減数選択法により、赤池の情報量基準（AIC）が最小となる変数の組合せを把握した。なお、相関係数の絶対値が 0.7 以上の変数は一方のみを使用して多重共線性の影響を避けた。さらに、得られた回帰式を用いて対象地全体の地上部バイオマスの分布を面的に推定した。
## 6 ）モデルの妥当性の検討

本項では、統計解析に用いたプロット数が40数ヶ所と少なく、独立したデータセットを用いた検証が難しい。そこで、CARTについてはトレーニングデータの何割が正しく分類されたかを用い、一般化線形モデルについては決定係数を用いてモデルの妥当性を評価した。なお、inverseとlogのリンク関数を用いた場合は決定係数が求められない。そこで、得られた回帰式から求めた地上部バイオマスの推定値と実測値（アロメトリ式から求めた地上部バイオマスの算出値）を相関係数や散布図を用いて比較することで代用した。また、森林管理の面からは、地上部バイオマスの値を正確に推定するよりも、地上部バイオマスの大きさをいくつかのクラスに分けて森林劣化の程度を段階的に示した地図があると有用である。そこで、一般化線形モデルの結果について、地上部バイオマスを200 t/ha毎にクラス分けした場合の実測値と推定値の対応関係も調べた。さらに、対象地全体で推定した地上部バイオマスの分布図の妥当性を調べるために、サバ州林業局の森林管理区毎の報告書 $ ^{6,7} $ に記載された、林班毎の情報が記載された地図（保護林や劣化林の分布図）と比較し、推定結果の妥当性を確認した。

## （７）時系列高分解能衛星データによる地上部バイオマス変化のモニタリング

時系列での地上部バイオマスの変化を広域でモニタリングするには、リモートセンシングデータが有効であるが、森林の施業履歴などによる劣化をとらえる高分解能衛星データの過去のアーカイブは非常に限られる。これは、高分解能衛星は商用衛星であるため、注文による撮影以外のデータが極端に少ないためである。このため過去と現在の比較を行うのに、異なる衛星のデータを使用しなければならないことが多い。この場合、衛星センサの地上分解能のほか、観測波長帯が異なること（表(1)-1）にも留意が必要である。とりわけオブジェクトベース画像分類を用いる場合、地上分解能の違いに応じた分類パラメータの選択が必要となる。

ここでは地上部バイオマスの分布を調べた2009年のIKONOS衛星データと2003年のSPOT衛星5号のデータを用いて、地上部バイオマスの変化をモニタリングした。まず、バイオマスの変化をモニタリングするためには、2時期の衛星データの差分をとる必要があるため、IKONOS衛星データを参照して、SPOT衛星データの精密幾何補正を行った。さらにオブジェクトベース画像分類ソフトウェアのeCognition Developerを用いて、林相の類似したエリアの領域分割を行った。その際、時系列データの地上分解能が異なることから、それぞれのデータ間で平均オブジェクトサイズがほぼ同じサイズになるよう領域分割パラメータを調整した。次に、先のIKONOS衛星データ解析で得られたバイオマス推定を教師データとして過去のSPOT衛星データに対してバイオマスに関する画像分類を行い、過去（2003年）のバイオマス分布図を作成した。これを比較することにより2時期の地上部バイオマス変化を推定した。

## 4 ．結果及び考察

## （１）地上調査による樹高と地上バイオマスとの関係

プロットごとに計算した平均樹高、標準偏差、最大樹高、20mを超える木の本数と地上部バイオマスの関係を図示した。これらの中で地上部バイオマスと最大樹高(r=0.64)、そして20m以上の樹木の本数(r=0.77)には相関が見られた（図(1)-13、表(1)-5）。平均やメディアンは関連が見られ
るもののレンジの幅が小さかった。これら樹高から得られる変数を用いてクラスタ分析を行ったところ、得られたクラスタは地上部バイオマス量と対応していた（図(1)-14）。その中でバイオマス量が多いクラスタは最大樹高が高く、また最大樹高が大きなクラスタは、20mを超える樹木の本数によりクラスをさらに分けることが可能であった。このことから最大樹高と20mを超える樹木の本数はプロットを地上部バイオマスに応じたクラスに分けるのに有効であり、プロット内の最大樹高および20m以上の樹木の本数を求めることでプロットを地上部バイオマスに応じたクラスに分類することが可能であることが明らかとなった。このことは衛星LiDARから求められる最大樹高や林分の構造に応じた変数を用いることにより森林を地上部バイオマスに応じたクラスの分けることができることを示している。

表(1)-6 に地上部バイオマス量を目的変数、樹高関連変数を説明変数とした重回帰分析の結果を示す。地上部バイオマス量の推定には、樹高の中位数、平均値、樹高 20m 以上の本数が有意であることが示され、決定係数は 0.56 であった。各変数と地上部バイオマス量との相関は、中位数が -0.63、平均値が 0.55、樹高 20m 以上の本数が 0.42 であった。この結果より、樹高から得られる変数である中位数、平均値、樹高 20m 以上の樹木の本数は、直径から Brown の式で求められる地上部バイオマス量との相関が強く、この 3 つの変数からから地上部バイオマス量が推定可能であることが示された。中位数や平均樹高あるいは樹高 20m 以上の（いわゆる樹冠から突出している樹木）本数は、これに関係するパラメータを衛星 LiDAR や高解像度衛星画像から導出できる。

DBH-H 関係を指数モデルと非線形モデルにあてはめ結果を表(1)-7 に示す。表(1)-8 には、それぞれのあてはめ結果についての AIC を示す。図(1)-15 と図(1)-16 は指数モデルおよび非線形モデルにあてはめた結果のグラフである。指数モデルと非線形モデルの比較では、いずれの地上部バイオマス量のクラスにおいても、指数モデルと比較して非線形モデルの方が AIC の値が小さく、あてはめの結果が良かった。次に地上部バイオマス量クラス間を比較すると、指数モデルでは低バイオマス林と中バイオマス林では、定数項 a、係数 b ともにほぼ同じ値であり、高バイオマス林で、定数項 a が小さくなり、係数 b がやや大きくなった。すべての樹木を用いてあてはめた結果では、定数項 a、係数 b は、低・中バイオマス林と高バイオマス林の中間の値を示した。非線形モデルにおいても、低・中バイオマス林では、係数 a、定数項 b のいずれもほぼ同じ値であり、高バイオマス林において定数項 b が大きくなっていた。すべての樹木を用いてあてはめた結果では、やはり係数 a、定数項 b は、低・中バイオマス林と高バイオマス林の中間の値を示した。非線形モデルでは定数項 b は低・中バイオマス林で小さく、高バイオマス林では大きくなっており、すべての樹木を用いてあてはめた結果では中間的な値となっていた。モデルの選択としては、非線形モデルの方がややあてはまりが良かったが、それほど重要ではなかった。一方で全ての樹木を用いたモデルに比べて、バイオマをスクラスに分けてモデル化した場合、当てはまりがよくなった。これはバイオマスクラス、すなわち森林劣化に応じて D-H 関係が異なっており、それぞれに適当なモデルをあてはめる必要があることを示している。また低バイオマス林と中バイオマス林では係数の値の大きさに違いが少なかったので、実際にデータに適用する際には、低・中バイオマス林と高バイオマス林の2つのクラスに分けて予測を行うのが適当である。
図(1)-13 地上部バイオマスと樹高（H）からの変数の関係

図(1)-14 樹高から得られる変数を用いたクラスタ分析結果
図(1)-15 指数モデルへあてはめ結果

赤点、赤線：低バイオマス林、緑点、緑線：中バイオマス林、青点、青線：高バイオマス林、黒線：すべての樹木

図(1)-16 非線形モデルへのあてはめ結果
表(1)-5 バイオマスと樹高関連変数との関係

最大樹高	樹高中位数	樹高
標準偏差	20m以上の
本数	地上部
バイオマス
最大樹高	1.000	-	-	-	-
樹高中位数	0.147	1.000	-	-	-
樹高
標準偏差	0.616	0.002	1.000	-	-
20m以上の
本数	0.376	0.650	0.075	1.000	-
地上部
バイオマス	0.639	0.428	0.153	0.772	1.000

表(1)-6 地上部バイオマスを目的変数、樹高を説明変数とした回帰分析結果

(偏回幅係数・信頼区間等)

偏回帰係数の有意性の検定
変 数	偏回帰係数	標準誤差	標準偏回帰係数	F 值	t 值	P 值	判 定
中位数	-87.4486	24.3045	-1.1078	12.9459	-3.5980	0.0018	**
平均值	75.7241	25.7177	0.9540	8.6697	2.9444	0.0080	**
20m以上の本数	12.2381	5.8584	0.6052	4.3638	2.0890	0.0497	*
定数項	110.6979	317.0894	0.1219	0.3491	0.7307

表(1)-7 D-H 関係のモデルへのあてはめ結果

Model	対象プロット	係数	推定値	標準誤差	t value	Pr(>|t|)
Exponential	低バイオマス林	a	-15.7719	1.1684	-13.50	<2e-16 ***
b	11.3152	0.4006	28.25	<2e-16 ***
Exponential	中バイオマス林	a	-15.2990	1.2950	-11.82	<2e-16 ***
b	11.3530	0.4410	25.74	<2e-16 ***
Exponential	高バイオマス林	a	-19.6922	1.3344	-14.76	<2e-16 ***
b	12.8456	0.4424	29.03	<2e-16 ***
Exponential	すべてのプロット	a	-16.8695	0.7205	-23.41	<2e-16 ***
b	11.8291	0.2440	48.48	<2e-16 ***
Nonlinear	低バイオマス林	a	1.3257	0.0458	28.98	<2e-16 ***
b	53.4005	2.8786	18.55	<2e-16 ***
Nonlinear	中バイオマス林	a	1.4110	0.0507	27.81	<2e-16 ***
b	52.5126	2.6777	19.61	<2e-16 ***
Nonlinear	高バイオマス林	a	1.3204	0.0473	27.91	<2e-16 ***
b	58.6107	2.8447	20.60	<2e-16 ***
Nonlinear	すべての樹木	a	1.3522	0.0268	50.38	<2e-16 ***
b	54.8413	1.5760	34.80	<2e-16 ***