人工衛星データを利用した我が国の森林炭素蓄積量を推定するための機械学習モデル開発
代表研究者

加 藤 知 道北海道大学 大学院農学研究院 准教授

1 はじめに森林が炭素（= CO2）を、いつ・どこで・どのくらい蓄積しているかを調べることは、地球の気候変動の正確な予測に貢献し（IPCC, 2013）、持続可能社会を形成する上で非常に重要である。特に国土の 70%を森林で覆われている我が国の炭素蓄積量は、第二次世界大戦後の拡大造林と木材利用の減少により、森林バイオマスが増加しつづけているため、現在は大きな炭素吸収源であると考えられている（Fang, Kato et al., 2014,
PNAS など）。そのため、京都議定書では、日本の排出量削減目標 6%のうち 3.8%を森林に割り当てられているなど、大きな生態系サービスを提供するものと非常に期待されているが、その詳細な時間・空間分布はわかっていない。
それらを克服するために、これまで、人工衛星データを利用した森林炭素蓄積量の推定が、世界の各地域で盛んに試みられている（Saatchi et al., 2008, Remote Sens. Env）
。特に、機械学習モデルを利用した方法が発達し始めているが、一方で学習用の地上検証データを入手することが非常に困難なために、研究が実現した地域が限られている。
そこで著者らは、森林地上部構造を強く反映する合成開口レーダー（PLASAR: JAXA）等の衛星データを入力値とした機械学習モデルを構築し、我が国の森林地上部炭素蓄積量のマップを作成することを目指している。我々はその初めの一歩として、高解像度の航空機光検出と測距(Lidar)データに基づき計算された森林地上部バイオマス量を学習データとし、茨城県のスギ・ヒノキ林および大分県の落葉広葉樹林を対象として、地域スケールの森林地上部バイオマス量の推定を行うランダムフォレスト機械学習モデルの開発を目的とした。

2 森林地上部炭素蓄積量の推定
2-1 茨城県のスギ・ヒノキ林における森林地上部炭素蓄積量の推定（１）はじめに茨城県北部（図１）において、日本国土の 18％と 10%を占める主要な植林品種であるスギ(Japanese cedar)とヒノキ(Japanese cypress)を対象として森林地上部炭素蓄積量（地上部バイオマス）の推定を行った。

図 1. 研究対象（茨城県北部のスギ(Japanese cedar)とヒノキ(Japanese cypress)）1
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

（２）データの整備
データ解析の流れとしては、a) 衛星データ収集と加工（リサンプリング、統一した座標系の適用、フィルタリングと画像くり抜き）、b) PALSAR-2 (JAXA)衛星データより景観テクスチャの抜き出し、c) 各種衛星画像からの複数指標の計算（HV & HH, NDVI, EVI）d) 機械学習モデルの構築（最適変数の選別、ハイパーパラメタの調整）、e) 茨城県の自治体における検証となる（図２）。
機械学習モデルは、航空機 Lidar(レーザー)観測に基づいた森林地上部バイオマス（=森林地上部炭素蓄積量）データを学習データとして利用する。本研究では茨城県庁が測量会社に発注した当該航空機観測（2020年 7 月実施）のデータを、各種手続きを経て入手し利用した。測量会社は、Lidar による波形解析で、地面高(DTM: Digital Terrain Model)と表面高(DSM: Digital Surface Model)の間の差分として、群落高(DCHM:Digital Canopy Height Model)を計算し、樹木個体の頂点での高さを樹高として、あらかじめ地上計測で求められたアロメトリー式（樹高や胸高直径 DBH からバイオマス量を求める経験式）を使用して、個体の材積量（材木生産の対象となる幹の容量）を計算している。対象個体は、スギで 201,854 本、ヒノキで 69,374本になった（表１）。これら個体材積量に材積密度(WD: Wood density)と拡大係数（BEF: Biomass ExpansionFactor）をかけることで地上部バイオマス量（=地上部炭素蓄積量）に変換することができる。これら処理を施し、20m x 20m メッシュサイズの地上部バイオマス量の学習データを得ることができた。

図 2. データ解析のフロー. AGB, 地上部バイオマス(=地上部炭素蓄積量); GLCM, gray-level co-occurrence
matrix; HH, horizontal transmit–horizontal channel; HV, horizontal transmit–vertical channel; MSI,multispectral instrument; VH, vertical transmit–horizontal channel; VV, vertical transmit–verticalchannel.

2
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

表 1. 航空機 Lidar から得られた樹木計測に関する統計値
SpeciesStand variable

Mean

Japanese cedar

Tree height (m)Diameter at breast height(cm)Stem volume (m3 ha−1)
Biomass (Mg ha−1)
Tree height (m)Diameter at breast height(cm)Stem volume (m3 ha−1)
Biomass (Mg ha−1)

Japanese cypress

Min.

Max.

24.1
24.1

Standarddeviation
5.2
5.32.19.946.7
78.0

403.6
155.9
19.2
27.0

170.7
65.9
4.3
5.8

0.3
0.1
2.3
10.1

1516.3
585.6
39.6
72.0

585.9
295.7

246.2
124.3

0.3
0.1

1800.0908.4

Samplesize
201 854

69 374

機械学習モデルの入力には、森林地上部構造をよく反映すると考えられている日本の宇宙航空研究開発機構 JAXA による Advanced Land Observing Satellite-2 (ALOS-2) 衛星に搭載の合成開口レーダ Phased Array
L-band Synthetic Aperture Radar-2 (PALSAR-2)を利用した。4 つの偏光モード(HH, HV, VH, VV)における後方散乱係数g0 を以下のように求めた。

（１）
ここで DN はデジタルナンバー、CF はキャリブレーション係数で-83 になっている。さらに 3x3 ピクセルを対象とした LEE スペックルフィルターを施し、林分上の不均一性によるスペックルノイズを減少させた。その後 3x3 ウィンドウサイズの Gray-level co-occurrence matrix (GLCM)（相対修正ベクター: d = 1, θ = 45°）における HV と VH のレーダー係数及びテクスチャ（各種統計値：表２）を計算した。

表２. PALSAR データのリスト
Variables (Abbreviation)
HVHH
Polarization
VV
VH
I1 [50]I2 [51]I3 [52]I4 [53]
Radar IndicesI5 [50]I6 [51]I7 [52]I8 [53]I9 [54]
Texture (HV, VH)

Definition
Horizontal transmit-vertical channel
Horizontal transmit-horizontal channel
Vertical transmit-vertical channel
Vertical transmit-horizontal channel
HH-HV
HV+HH
(HH-HV)/(HV+HH)
HV/HH
VH-VV
VH+VV
(VH-VV)/(VH+VV)
VH/VV
8*HV/(HH+VV+2*HV)Mean (ME)3
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

Variance (VA)
Homogeneity (HO)Contrast (CON)Dissimilarity (DIS)
Entropy (ENT)
Second Moment (SM)
Correlation (COR)

さらに欧州宇宙機関(ESA)の衛星である Sentinel2 による光学センサーのデータも補助として利用した。
我々は 2019 年 5 月 8 日の Level2-A の Band2-12 までのデータをダウンロードし、３x３ピクセルの平均を行い、表 3 に示す植生指標(Vegetation Indices)を計算した。NDVI, EVI 等の植生指標は、地上植生の葉の総量や光合成活性を表すものとして、これまでたくさんの研究によって利用されてきており、森林構造を反映する PALSAR データを補完するものとして利用した。

表 3. 入力に利用した Sentine2-MSI データおよび計算した植生指標
Variables
Bands, indices (Abbreviation)

Multispectral Bands

Definition (Central Wavelength)

Band2(B2)
Band3(B3)
Band4(B4)
Band5(B5)
Band6(B6)
Band7(B7)
Band8(B8)
Band8A(B8a)
Band11(B11)
Band12(B12)Blue,490nm
Green,560nm
Red,665nm
Red edge,705nm
Red edge,749nm
Red edge,783nm
Near Infrared (NIR), 842nm
Near Infrared (NIR), 865nm
SWIR-1, 1610nm
SWIR-2, 2190nm

NDVI [55]
EVI [56]
Vegetation Indices

DVI [57]
GARI [58]
SAVI [59]4公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

GNDVI [60]
GDVI [61]SR [62]GRVI [63]（３）機械学習モデル
本研究ではランダムフォレスト(Random Forest)を利用した。コンピュータ言語 Python で sklearn パッケージを使用し、48 種類のデータ（10: Sentinel2-MSI スペクトルバンド, 9: Sentinel-2 MSI 植生指標, 4:PALSAR2 後方散乱係数, 9 PALSAR2 レーダー後方散乱係数指標, 16: PALSAR2 テクスチャ）を入力とし、航空機ライダー地上部バイオマスデータを 80%, 10%, 10%に分割し、それぞれ学習、検証、テストサンプリングに使用した。ランダムフォレストは、たくさんの決定木を構築し、変数ごとの重要性を確認しながら、学習データとしてのバイオマスデータを再現する最良のモデルを求める。モデルの評価には、平均平方誤差 RMSE、R2、平均絶対誤差 MAE、相対 RMSE を利用した。学習データの中で、250 Mg/ha 以上の非常に大きなバイオマス量を測定エラーとして除去した。

（４）結果と考察
ランダムフォレストモデルでは、入力変数のカテゴリー数と、説明変数の量に影響を受けて、入力変数の重要度（importance）が変化するが、入力変数の量を増やすと計算時間がかかりすぎてしまうという問題がある。そこで我々は PALSAR-2 と Sentinel-2 の 2 つのカテゴリーごとに Gini 係数が 0.05 以上の変数のみを利用して、計算を行うことにした。
その結果は、
スギの PALSAR-2 については、
VH mean、VH variance、HV variance、VH correlation、HV correlation が、スギの Sentinel-2 については Band 5、Band 9、Band 8a、Band 11、SR、NDVI、Band 12 が、ヒノキの PALSAR-2 については、VH mean、VH variance、HV mean、HV variance が、ヒノキの Sentinel-2 については Band 12、Band 4、Band 9、Band 11、Band 5、Band 8a、Band 6 が重要な変数として選択された。
ランダムフォレストアルゴリズムにおいて、ハイパーパラメータの調整は、モデルの性能を決める重要な課題である。繰り返し計算の結果、Maximum depth of the decision tree (MD)、Minimum number of samples(MS) 、Minimum number of samples required to be at every leaf node (ML)、Number of trees in the forest(EST)を、それぞれ 10,2,1,200 に設定した。
構築したモデルは、スギ林に対して作成したモデルは、R2 = 0.31, root-mean-square error (RMSE) = 54.38
Mg ha−1, mean absolute error (MAE) = 40.98 Mg ha−1, relative RMSE (rRMSE)= 0.35, ヒノキ林に対して作成したモデルは R2 = 0.37, RMSE = 98.63 Mg ha−1, MAE = 76.97 Mg ha−1, rRMSE = 0.33 のある程度高い精度を示した（図４）。それらモデルを最終的に茨城県全体に対して適用した結果、17 の県内市町村におけるスギの地上部バイオマス量は 5.27 Pg で、単位面積あたりの平均値は 146.50 Mg ha−1 、標準偏差は 44.37
Mg ha−1、ヒノキ林は、3.56 Pg、293.12 Mg ha−1 、78.48 Mg ha−1 であった（図５）。
我々が計算した地上部バイオマス量を、従来からの方法である森林簿（対象林小班における樹種・林齢を変数とする経験式）から求めたものについて、17 の市町村（水戸、日立、石岡、常陸丘、高萩、北茨城、笠間、ひたちなか、常陸宮、那珂、鉾田、小美玉、茨城、大洗、下里、東海、大子）の平均値を求めて比較した（図６）。その結果は、両者の相関係数は 0.99 と非常に高く、森林簿はスギで 25%の過大評価、ヒノキで53%の過小評価をしていることがわかった。
このように、航空機ライダーによる森林バイオマス量を学習データとし、衛星データを入力データとした
機械学習モデルは、広域で高い精度で我が国の主要な商業用針葉樹であるスギとヒノキの地上部バイオマス量を再現することができた。また森林簿による従来の推定方法は、正確ではない可能性が高いことも示唆した。

5
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

図 4. 茨城県における観測と衛星データによる推定された森林バイオマス（(a)スギ、(b)ヒノキ）

図 5. 茨城県における衛星データから推定された森林バイオマス（(a)スギ、(b)ヒノキ）

6
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

図 6. 茨城県北部１７市町村における衛星データと森林簿から推定された森林地上部バイオマス（(a)スギ、(b)ヒノキ）

2-2 大分県の落葉広葉樹林における森林炭素蓄積量の推定（１）はじめに
我が国の森林タイプの代表的なものとして落葉広葉樹林があるが衛星データを利用した森林地上部バイオマス（=炭素蓄積量）の推定はほとんど行われていない。そこで、大分県のクヌギに代表される落葉広葉樹林を対象（図７）とし、航空機ライダー観測に基づいた地上部バイオマス量を再現する機械学習モデルの構築を目的とした。

図７. 大分県における対象地（落葉広葉樹 DBF: 赤部分）

7
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

（２）データの整備入力に利用する衛星データとして用意した Sentinel-2(ESA)による各種植生指標(NDVI, EVI, DVI, SR)及び PALSAR-2(JAXA) L-band SAR 画像の他に、PALSAR-2-FNF(森林分布図: JAXA; Shimada et al., 2014), SRTM（スペースシャトルによる地形情報, NASA）、Sentinel-1(ESA)なども利用することにしたが、それらのデータ容量が非常に大きいためランダムフォレスト回帰モデルを操作する上での大きな妨げとなる。そのため本研究では主成分分析(PCA)を利用して複数種類のデータをいくつかの軸にまとめる作業を行い、データ処理時間を減少させることにした。3x3 平均フィルターののちに、5x5 の Lee フィルターによってスペックルノイズを低減させた。
航空機ライダー観測に基づいた落葉広葉樹林の森林地上部バイオマスを学習データに利用した。それらはライダーデータより求められた群落高(DCHM)に対して地上１１ヶ所で計測された値より構築したアロメトリー式を利用して求めた材積量に、材積密度 0.646 Mg/m3 及び拡大係数 1.33 をかけたものであり、10mx10m のメッシュで 89,324 ピクセルのデータ容量であった。地上部バイオマス密度は、0.1 から 701.5 Mg/ha の範囲をとり、平均 121.267Mg/ha で標準偏差は 60.47Mg/ha であった。

（３）機械学習モデル
ここでは入力データの種類の異なる 4 つのランダムフォレスト(RF)回帰モデルを試した。Model 1 –土地基礎情報 (PALSAR-2-FNF, SRTM, Sentinel-1) と PCA グループ（Sentinel-2 の PCA と植生指標のテクスチャ）,
Model 2 - 土地基礎情報と VI グループ(Sentinel-2 の植生指標とそのテクスチャ), Model 3 - 土地基礎情報とオリジナルバンドグループ(Sentinel-2 のマルチバンドデータとそのテクスチャ), All data -上記全て（表３）.
表３. ランダムフォレスト回帰モデルの概要

RF-based model
Model1Model2Model3All data

入力データ
土地基礎情報, 3 PCA bands, and 30 relative textures.
土地基礎情報, 10 VIs, and 30 relative textures.
土地基礎情報, 10 original bands, and 30 relative textures.
土地基礎情報、3 PCA bands, 10 VIs, and, 10 original bands, 90
relative textures.

入力変数の総数

707777150（４）結果と考察
4 つのランダムフォレストモデルを試したところ、Model1 が最も良いパフォーマンスを示した（図 8; R2 =0.69, root-mean-square error (RMSE) = 33.59 Mg/ha, mean absolute error (MAE) = 25.01 Mg/ha, relative
RMSE (rRMSE) = 0.27）。Recursive Feature Elimination (REF)プロセスによって選択された変数は、34 個のうち 19 個が PCA とそのテクスチャに関するものであった。
Model1 により 34 個に絞り込まれた変数でモデル構築し、計算したところ、落葉広葉樹林の森林地上部バイオマスは 32.69 から 437.72 Mg/ha の範囲をとり、平均 128.84Mg/ha で標準偏差は 32.24Mg/ha であった。
大分県での総計は 14,507,596.8 Mg に登った（図 9）。
本研究では、新たに PCA を利用して変数容量を減らすことなどにより効率的なモデルを構築するとともに、これまであまり挑戦されてこなかった、我が国での落葉広葉樹林の森林地上部バイオマス推定について成功した。

8
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

図 8: 大分県の落葉広葉樹林における地上部バイオマスの観測とモデルの比較: (a) model1, (b) model2, (c)model3, (d) All model.

図 9: 大 分 県 の 落 葉 広 葉 樹 林 に お け る 地 上 部 バ イ オ マ ス . (a) AGB distributions derived fromsatellite+RF-RFE (left) and airborne LiDAR (right), (b) histogram of mapping AGB in 10m grid cell,(c) AGB map of DBF in study area with close-up images (d) and (e).

9公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度

3 おわりに本研究では、茨城県のスギ・ヒノキ林及び大分県の落葉広葉樹林を対象とし、人工衛星を利用した森林地上部炭素蓄積量（＝森林地上部バイオマス）を推定するための機械学習モデルの構築に成功した。これまで我が国ではあまり行われてこなかった地域スケールでの推定であり、PCA によるデータの圧縮を利用した計算の効率化など手法的な新規性も見出すことができたことは非常に大きな成果であった。今後はより広域なスケールへ拡張していくべくすでに研究を進めているが、一方による広域化で樹木タイプが増えることによる情報要求量の増加や、その大容量データの効率的な処理について、さらなる技術的な開発必要であることが予想される。我が国は南北に長いための樹種タイプの豊富さや、地形の複雑さを克服して森林地上部炭素蓄積量の推定を目指すことは、森林に国土の 7 割が覆われる日本だけでなく同様な各国にとって非常に意義深いものであると考えられる。

【参考文献】
Saatchi S, Buermann W, ter Steege H et al (2008) Modeling distribution of Amazonian tree species
anddiversity
using
remotesensingmeasurements.
RemoteSens
Environ112(5):2000–2017. https://doi.org/10.1016/j.rse.2008.01.008
IPCC, 2013: Climate Change 2013: The Physical Science Basis.
Fang, J.Y., Kato, T., Guo, Z., Yang, YH., Hu, HF., Shen, HH, Zhao, X., Kishimoto-Mo, AW., Tang, YH.
Houghton, RA., 2014: Evidence for environmentally enhanced forest growth. Proc. Nat. Aca. Sci.,
111(26), 9527-9532, doi: 10.1073/pnas.1402333111
Shimada, M., Itoh, T., Motooka, T., Watanabe, M., Shiraishi, T., Thapa, R., Lucas, R., 2014. New
global forest/non-forest maps from ALOS PALSAR data (2007–2010). Remote Sensing ofEnvironment 155, 13-31.

〈発題名Using multisource satellite products toestimate forest aboveground biomass in
Oita prefecture: An approach withimproved accuracy and computational
efficiency
Estimation of Forest Aboveground
Biomass of Two Major Conifers in
Ibaraki Prefecture, Japan, fromPALSAR-2 and Sentinel-2 Data

表資料〉

掲載誌・学会名等

発表年月

Geocarto International

2022 年 10 月

Remote Sensing

2022 年 1 月

10
公益財団法人電気通信普及財団
研究調査助成報告書 第38号 2023年度