特集「環境リスクと統計解析―データ基盤構築と解析―」［研究ノート］

統計数理（2013）
第 61 巻 第 2 号 307–322
c 2013 統計数理研究所


人工林長期継続調査データを利用した林分成長モデルのパラメータ推定12

光田 靖 ・細田 和男 ・家原 敏郎

2（受付 2012 年 12 月 28 日；改訂 2013 年 11 月 5 日；採択 11 月 11 日）

要

旨

林業の収穫量を把握することを目的として長期間継続調査がなされている試験地（収穫試験地）のデータを用いて，光合成や呼吸による物質収支に基づくスギ人工林林分成長モデルのパラメータ推定を行った．本研究における林分成長モデルは 1）葉群による日射の吸収プロセス，2）
吸収した日射エネルギーを用いた光合成による生産プロセス， 3）光合成による日射エネルギー
変換効率が気温および湿度から律速される律速プロセス， 4）樹木個体を維持するための呼吸および 5）器官の付け替え（ターンオーバー）プロセス，および 6）余剰光合成生産物を各器官へ配分してサイズ成長を行う配分プロセスといった 6 つのプロセスから構成される．収穫試験地において観測された樹木サイズ（直径および樹高）から部位別（葉，幹，枝および根）バイオマスを推定し，さらに連続する調査データの差分としてバイオマス成長量をもとめた．また，試験地における調査期間の日射量，平均気温および湿度といった気象値を推定した．これらの部位別
バイオマスおよび気象値の時系列データを入力値とし，バイオマス成長量を目的変数として，
ベイジアンキャリブレーションによりパラメータ推定を行った．推定されたパラメータを用いてモデルの再現性を検証したところ，大きな過大推定をするケースもあったが，おおむね良好に林分成長量を再現できていた．
キーワード： パターン試行モデリング，ベイジアンキャリブレーション，スギ人工林．
1. はじめに林業という産業は木材を供給するという森林の生態系サービスを利用する産業であり，そのため自然条件に大きく左右されるという特徴がある．林業経営において環境条件に対する判断を誤ると，植栽した樹木が成長しない，土砂災害，風害や虫害といった自然災害が発生するといった自然環境条件に規定されるデメリットが生じる．将来的な気候変動はこれら林業経営におけるデメリット発生に大きく影響する可能性が高く，木材生産という生態系サービスに対する環境リスクとなる（松本 他, 2006; 松本, 2008; 農林水産技術会議事務局, 2011）．さらに，戦後の木材供給不足を解消するため広葉樹の林を伐採してスギやヒノキなど針葉樹の林へと転換
する拡大造林政策がとられ，森林の約 40％は針葉樹人工林となっている．これらの針葉樹人工林は樹種が単一（スギやヒノキ）であり，また構造的にも単純（樹齢が等しく，サイズも変動が1

森林総合研究所 四国支所：〒780–8077 高知県高知市朝倉西町 2-915 （現 宮崎大学 農学部：〒889–2155

宮崎県宮崎市学園木花台西 1 丁目 1）
2

森林総合研究所 森林管理研究領域：〒305–8687 茨城県つくば市松の里 1

308

統計数理

第 61 巻

第 2 号 2013

小さい）であるため，気候変動の影響が顕著になる可能性がある．自然環境条件に応じて生じる林業におけるデメリットは，間伐など森林の取り扱いで低減できる部分もあるが，大部分は自然環境条件によって規定されるため，短期的に対応できることは少なく，長期的視野に立った対応が求められる．気候変動という環境リスクに長期的に対応するためには，環境に応じた樹木の反応に関するデータベースおよびそれに基づいた環境リスク評価・予測手法の開発が必要不可欠である．
森林における環境リスクを評価するという目的ではないが，林業収益性の観点から立地環境に応じた樹木の成長を直接的な観測によって評価するため，長期的な継続調査試験地として「収穫試験地」が日本全国の国有林に設置されている（森林総合研究所, 1996）．収穫試験地とは国有林と独立行政法人森林総合研究所が協同で維持，管理している試験地であり，試験区内の樹木個体をナンバリングして個体識別しながら胸高直径（地面から 1.2 m もしくは 1.3 m の高さにおける幹の直径）や樹高を定期的（基本的に 5 年もしくは 10 年間隔）に計測している．観測項目は限定的であるものの，全国に 200 箇所以上と多数分散して存在しており，計測期間は長いもの
で 70 年以上にわたり，どのような地域で，どのように林が成長してきたのかを実証する貴重なデータとなっている．収穫試験地はどの程度の木材をどの程度の期間で生産することができるのかを調べるために設置されたものであるが，その継続調査データは様々な研究に応用することが可能である．例えば最近では，間伐を行った場合と行わなかった場合について森林の二酸化炭素吸収量にどのような差が生じるのかを明らかにする研究に利用されている（細田 他,
2012）．樹木の成長は過去の気象や個体間競争といった環境条件に対する応答の累積であり，例
えば過去の気象データとつき合わせることによって，収穫試験地における継続調査データは環境リスク研究における貴重な情報基盤となりうるだろう．
本研究においては，木材生産という生態系サービスに対する気候変動という環境リスクを予測するためのツールとして，森林の成長過程を気象環境への応答として記述する林分成長モデ
ルを，収穫試験地における長期観測データを用いて開発することを目的とする．なお，対象樹
種は日本の人工林において最も植栽面積の大きいスギとする．また，収穫試験地における継続調査データを用いて林分成長モデルのパラメータを推定することによって，収穫試験地データ
が環境リスク研究の情報基盤となる可能性について言及する．
2. 解析方法および使用データ
2.1 林分成長モデル
林分成長モデル研究には長い歴史があり，様々なモデルが提案されてきた．古くは林分蓄積と林齢の関係を曲線により回帰したものから，最近では樹木の生理的な反応の集合として成長を記述するモデルまで，そのコンセプトやスケールは多岐にわたる（Landsberg, 2003; Pretzsch
et al., 2008）．本研究においては，気象環境に対する樹木の反応を表現することを主題として，
光合成による物質生産と呼吸や器官の付け替えによる消費を対象とした物質収支を基盤とするモデル（物質収支モデル）を採用する．物質収支モデルも数多く提案されているが（例えば，Aberand Federer, 1992; Bergh et al., 2003），シンプルであるが現実的に適用しやすいモデルとして定評のある 3-PG モデル（Landsberg and Waring, 1997）をモデルの基本形として用いた．3-PG
モデルは日射エネルギー変換効率を用いて計算される総一次生産量を基礎として，気象条件から光合成を律速する物質収支モデルで，様々な樹種について適用事例がある（例えば，Coops etal., 2005）．
本研究においては実用性の観点から 3-PG モデルを簡略化および改良し，1）葉群による日射の吸収プロセス，2）吸収した日射エネルギーを用いた光合成による生産プロセス， 3）光合成による日射エネルギー変換効率が気温および湿度から律速される律速プロセス， 4）樹木個体を維人工林長期継続調査データを利用した林分成長モデルのパラメータ推定

図 1.

309

林分成長モデルの概要．

持するための呼吸および 5）器官の付け替え（ターンオーバー）プロセス，および 6）余剰光合成生産物を各器官へ配分してサイズ成長を行う配分プロセスといった 6 つのプロセスを持つモデルとした（図 1）．それぞれのプロセスは対象とする林分のサイズおよび気象要因によってその反応が決まるように設計されている．なお，本モデルは物質収支に基づくものであるため，林分のサイズを示すものとして各器官（葉，枝，幹および根）のバイオマス（単位面積あたりの乾燥重量 [Mg/ha]）を用いる．以下に各プロセスについて定式化を行う．
1）日射吸収プロセス
林分に入射した日射エネルギーは葉群により吸収され，光合成に利用される．この日射吸収
プロセスおよびそれに続く光合成生産プロセスのモデル化には門司 佐伯による葉群光合成モ
デル（Monsi and Saeki, 1953）の原理を適用した．本研究においては，5 層の葉群階層を設定して上層から順に各葉群階層における日射吸収量を計算し，入射した日射の減衰過程を表現する．
(2.1)

P ARi = P ARi−1 exp(−K × Bfi )(2.2)

AP ARi = P ARi−1 − P ARi

ここで，葉群階層 i を通過した後の光合成有効放射量（Photosynthetically Active Radiation: PAR[MJ/ha/month]）を P ARi とし，Bfi は葉群階層 i のヘクタールあたり葉量を示す．また，K は
門司 佐伯モデルにおける吸光係数であり，本研究においては樹種で一定とした．葉群階層 i において吸収された光合成有効放射量（Absorbed PAR, AP ARi ）は，入射した放射量と通過後の放射量との差分として計算される．光合成有効放射量は入力気象値である全天日射量の半分とし，林分に入射する前の光合成有効放射量を P AR0 とする．本来の門司 佐伯モデルにおいては葉面積指数を用いるが，本モデルにおいては各器官のバイオマスを基礎単位としてモデル計
算を進めるため便宜的に葉量を用いた．
葉群を単相で取り扱う 3-PG モデルに対して，本モデルにおいては葉群の階層構造を新たに導入した．これは 3-PG モデルが総一次生産量の計算に固定の日射変換係数を用いるのに対し，本モデルでは光強度に応じた光合成速度によって総一次生産量を計算することに対応している．
これらの改変は，光強度の高い林冠上部では光合成速度が高くなり，光強度の低い林冠下部では光合成速度が低くなるという，現実に観測される現象を表現するためのものである．このことによって，個体サイズ分布の違いによって林冠構造に違いが生じ，光合成生産量が変わるこ

310

統計数理

第 61 巻

第 2 号 2013とを表現できるようになる．
2）光合成生産プロセス
各葉群において吸収した日射エネルギーによって光合成による物質生産が行われる．門司 佐伯による葉群光合成モデルにおいては，光合成速度は入射する光強度によって変化する．先述の日射吸収プロセスによって表現されるように各葉群階層に入射する日射は減衰することから，
光強度と光合成速度の関係をモデル化して各葉群階層における光合成速度を推定する．ここでは，モデルの変数となるのが月間光合成有効放射量であることから，月間光合成有効放射量を
光強度，月間の放射量から光合成生産量への変換効率を光合成速度とみなして，月間放射量と光合成速度の関係を非直角双曲線によって次式のように記述し，光 光合成曲線とする．

φP ARi−1 + Pmax − (φP ARi−1 + Pmax )2 − 4φP ARi−1 θPmax
(2.3)
Pi =
2θ
ここで，Pi は葉群階層 i における光合成速度（[Mg/MJ]) を示し，Pmax は最大光合成速度を示
すパラメータ，φ および θ は曲線形を規定するパラメータである．なお，Pmax および φ は非負の値をとり，θ は 0 から 1 の値をとる．葉群階層 i におけるストレスのない状態での光合成によるポテンシャル物質生産量（P Pi [Mg/ha/month]）は日射吸収プロセスにより計算される日
射吸収量と光合成速度から計算される．
(2.4)

P Pi = Pi × AP ARi3）光合成律速プロセス
光合成生産プロセスにおいて用いた光 光合成曲線は環境ストレスが全くない条件下での光合成速度を与えるものである．光合成速度は光強度だけでなく，様々な要因により律速されるが，本研究においては 3-PG モデルで用いられる律速因子のうち温度と湿度を律速因子として用いた（Landsberg and Waring, 1997）．光合成速度は低温で光合成酵素の活性が低下すること
で下がり，気温の上昇とともに上昇していくことから，0 から 1 の値をとる温度律速係数（MT ）を次の式で定義する．
1
(2.5)
MT =
1 + exp(−β1 (T − β2 ))ここで，T は平均気温（[◦ C]）を示し，β1 および β2 はパラメータであり，非負の値をとる．
また，大気の湿度が下がって乾燥することによって葉の気孔が閉じて光合成速度が低下する．
ここでは大気の乾燥を示す因子として相対湿度と気温から計算される大気飽差を用い，0 から1 の値をとる大気飽差律速係数（MV ）を次の式で定義する．
exp(−β3 (V P D − β4 )) (V P D > β4 )
(2.6)
MV =
1
ここで，V P D は大気飽差（Vapor Pressure Deficit: VPD[kPa]）を示し，β3 および β4 はパラメータであり，非負の値をとる．大気飽差の計算には，乾燥の指標として強調するため，気象入力
値である平均相対湿度と最高気温を用いた．
各葉群階層ごとに光 光合成曲線により算出されたポテンシャル光合成速度を気温および大気飽差律速係数によって補正して，気象要因によって律速された実効光合成速度と各葉群階層で吸収された日射量とをそれぞれ掛け合わせて各葉群階層における光合成生産量を算出する．
各葉群階層の光合成生産量を合計して林分の総一次生産量（Gross Primary Production: GPP）
を求める．
(2.7)

GP P =5
i=1

P Pi × MT × MV

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定3114）呼吸プロセス
光合成によって生じた生産物は各器官において呼吸により消費され，その量はそのバイオマ
スに比例する．呼吸量とバイオマスとの関係式はしばしば次式のようにべき乗式で表現される（Mori et al., 2010 参照）．
(2.8)

rs

Rj = rrj,T × Bj j

ここで，Rj は月間呼吸量（[Mg/ha/month]），Bj はバイオマス，rrj,T は月平均気温によって変動する呼吸速度パラメータ，および rsj はスケーリングパラメータを示し，添字 j は部位（葉群：f ，幹：s，枝：b，および根：r ）を示す．本研究においてはスケーリングパラメータを葉群については 1(rsf = 1）とし，幹，枝および根については 3/4(rss = rsb = rsr = 3/4）とした．
以前行った小サンプルでの予備的な解析では呼吸速度パラメータを一定としていた（Mitsuda
et al., 2010）．しかし，呼吸速度は気温により上昇することが知られており，気候変動の影響を
より現実的に推定するためには，呼吸速度と気温の関係をモデルに組み込む必要がある．そこ
で，次式により気温による呼吸速度の律速を導入した（千葉, 2004）．
(2.9)

rrj,T = β5 exp(β6 T )ここで，T は平均気温（[◦ C]）を示し， β5 および β6 はパラメータであり，非負の値をとる．呼吸速度推定式パラメータは葉群（β5 F , β6 F ）と木部（幹，枝および根， β5 C, β6 C ）に分けて推定する．
林分の総一次生産量から総呼吸量を差し引いて，林分の純一次生産量（Net Primary Production:
NPP）が次式により計算される．
(2.10)N P P = GP P −Rj
j∈{f,s,b,r}5）ターンオーバープロセス
スギやヒノキなどの常緑針葉樹においても，古い葉を落として新しい葉へ付け替えを行っている．同様に枝や根においても古い部位から新しい部位へのターンオーバーがある．なお，幹についても樹皮のターンオーバーがあるが，本研究においては対象外とした．ターンオーバー
量は各器官のバイオマスに比例し，季節性があるため月単位ではなく年単位で次式により計算される．
(2.11)

Lj = lrj × Bj

ここで，Lj は部位 j における年間ターンオーバー量（[Mg/ha/year]），lrj は年間ターンオーバー率を示す．根については細根（f r ）と粗大根（cr ）に分け（バイオマスの比率は 1:99 と仮定），そ
れぞれについてのターンオーバー率（lrf r および lrcr ）を設定する．
年間純一次生産量（月間 NPP の年間積算，N P PA ）から総年間ターンオーバー量を差し引いて，余剰光合成産物量（すなわち年間バイオマス成長量 [Mg/ha/year]）が次式により計算される．
(2.12)
G = N P PA −Ljj∈{f,b,f r,cr}6）余剰光合成産物配分プロセス
余剰生産物を各器官に配分して，各器官のサイズを成長させる．本研究においては樹木個体の総バイオマスに対する各器官バイオマスの比は， Fukuda ら（2003）の部位別バイオマス割合
推定モデルを利用して樹齢の関数により計算される．

統計数理

312
表 1.

第 61 巻

第 2 号 2013

収穫試験地の概要．調査区間数は解析に用いた区間数を示す．詳細は本文参照．

以上のように，モデルの入力となるのは各器官のバイオマスおよび気象因子として全天日
射量（月間積算値，光合成有効放射量へ換算してモデルへ入力），平均気温，最高気温および平均相対湿度（いずれも月平均値）である．モデル計算におけるタイムステップは日射吸収，光合
成生産，光合成律速および呼吸プロセスでは月単位とし，月単位の NPP を合計して年単位のNPP を計算した後，ターンオーバーおよび余剰光合成産物配分プロセスの計算を年単位で行う．
2.2 使用データ
モデルの入力値は林分の器官別バイオマスおよび月別気象値（全天日射量，平均気温，最高気温および平均相対湿度）であるため，モデルパラメータを推定するためにはこれらの実測値が必要となる．本研究においては，林分バイオマスの時系列データとして収穫試験地の継続調
査データを用い，気象値の時系列データとして 30 年間の平年値および近傍気象観測点の時系
列データを利用した．本研究で用いる収穫試験地は関東森林管理局管内に設定された 7 箇所のスギ人工林試験地で，間伐区と無間伐区が設定された試験地もあることから，合計で 13 試験区の時系列データを用いた．各試験区の諸元を表 1 に示す．モデルパラメータ推定には林分バ
イオマスの初期値と成長量が必要であることから，ある調査時点のデータとその次の調査時点のデータを 1 つの調査区間データとして取り扱い，調査区間期首のデータを林分バイオマスの初期値，期首と期末の林分バイオマスにおける差分を林分成長量として取り扱った．なお，最近傍気象観測点のデータ取得期間が足りない場合があり，気象データのない調査区間のデータは利用しなかった．また樹高を全数計測した場合とサンプル測定した場合があり，期首と期末で樹高測定方法が異なる調査区間も対象外として，合計で 61 調査期間のデータを解析に用いた．全調査区間データのうち，モデルのパラメータ推定に 39 データを用い，残りの 22 データ
はモデルの検証に用いた．
収穫試験地における測定項目は個体の直径および樹高であるため，バイオマスに変換する必要がある．樹高を全数調査していない場合については，各試験区の各調査時点それぞれにおいてサンプル調査された樹高と直径から関係式を作成し，未測定木の樹高を推定した．各個体について樹高と直径から細田ら（2010）の手法を用いて幹材積を推定した．次に京都議定書報告でも用いられているスギの容積密度数（0.314 [Mg/m3 ]）を用いて幹材積から幹バイオマスに変換した．調査時点の林齢から Fukuda ら（2003）の部位別バイオマス割合推定モデルを用いて部位

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定

313

別割合を計算し，幹バイオマスから各部位のバイオマスを計算した．各個体の部位別バイオマ
スを積算して林分の部位別バイオマスを求めた．また，各調査区間の期首と期末の差分からバイオマス成長量を計算した．
各試験区の各調査区間において期首の測定時から期末の測定時まで，月単位で各気象値の推定を行った．まず，各試験区の最近接気象観測点において，気象庁発行の数値メッシュ気候値
2000 から得られる 30 年間の平年値を基準として，調査区間に該当する時系列気象観測値について，基準値からの変動を正規分布の累積密度関数を利用して計算した．試験区が存在するメッ
シュの 30 年平年値を基準として，計算された変動より試験区における気象値の時系列データを推定した．なお，メッシュ気候値 2000 には相対湿度が含まれていないため，志賀（2003）の方法により水蒸気圧を介して相対湿度を推定した．
2.3 パラメータ推定
収穫試験地の継続調査データおよび時系列気象データを用いて，ベイジアンキャリブレー
ションによって林分成長モデルのパラメータを推定した．ベイジアンキャリブレーションとはマルコフ連鎖モンテカルロ（MCMC）によるベイズ推定によって，観測データからモデルのパラメータを推定する方法である（山道・角谷, 2009）．本研究において観測値は樹木の個体サイズのみであり，その後背にある光合成や呼吸といったプロセスについての観測値はない．先述のように，環境リスクを評価するためには環境に対する樹木の対応プロセスをモデルに組み込む必要があるが，現地観測を行なって各プロセスの反応をパラメータ化することは難しい．さらに，光合成などの生理反応を個葉で計測することは比較的容易であるが，個体単位さらに林
分単位で計測するためには大規模な実験設備が必要となる（例えば, 大谷, 2008）．このような観測の難しいプロセスについて，そのプロセスがもたらす帰結（ここでは林分のバイオマス成長量）からプロセスのパラメータ化を試みるためにはパターン指向モデリングの考え方が有効である（山道・角谷, 2009）．本モデルのようにパラメータ数が多くなると，総当たり的にパラメータを動かして組み合わせを試すキャリブレーション法は現実的ではなくなり，その実際的な解決手法として既往の知見を活用できるベイジアンキャリブレーションが有効である．
将来的な気候変動が林分成長に及ぼす環境リスクを予測するための林分成長モデルを開発することを目的としていることから，本研究においては特に気温への反応に着目し，光合成律速プロセスにおける気温律速係数および呼吸プロセスにおける呼吸速度の気温律速についてモ
デルのパラメータ推定を行った．本研究において取り扱う環境因子は気象要因（日射量，気温および湿度）のみであるが，現実的には土壌や地形といった要因が強く影響している（松本 他,
2006）．立地条件のバラツキを考えると，本研究に用いるデータは 7 箇所の試験地から得られたデータであり，仮に試験区を別立地条件として考えても 13 試験区しかないので，土壌要因
や地形要因の影響をモデルに組み込むのは難しい．さらには植栽されたスギの遺伝的形質における差違が，林分成長に影響を与えている可能性もある．そこで，モデルの入力値となる林分の初期サイズ，日射量，気温，および湿度以外の要因によって生じる林分成長のバラツキを，試験区の差として捉えることとした．具体的には，最も試験区の差を反映すると考えられる光光合成曲線において最大光合成速度を意味するパラメータ Pmax に階層構造を導入し，スギの樹種平均的なパラメータ（Pmax.sp ）から各試験区固有のパラメータ（Pmax.plot ）を派生させる
階層ベイズモデルとした．このような階層化を行う利点は，スギの樹種平均的なパラメータである Pmax.sp をより適正に推定できる点にある．ここでは試験区差を考慮していない他のパラメータも含めて，スギという樹種の代表的な振る舞いを記述するようなモデルのパラメータ推定に集中することとする．本研究においては 13 試験区のデータを用いているが，立地条件が極端に良いあるいは悪い試験区がある場合，階層化を行わないとスギの樹種平均的なパラメー

統計数理

314
表 2.

第 61 巻

第 2 号 2013

林分成長モデルにおけるパラメータの概要．ただし，sd =バーを意味する．



1/tau．T.O. はターンオー

タの推定に偏りを生じさせてしまう．さらに，実際には直径と樹高の観測から推定する葉量だけでは捉えきれない林冠構造の時系列変化など，情報のない要因によって同じ試験区でも光合成速度が調査区間によって変化する可能性もあることから，よりよいパラメータ推定のために調査区間固有のパラメータ（Pmax.prd ）を導入した．試験区 i および試験区 i における調査区間
j についての固有パラメータは，上位のパラメータを平均値とする正規分布を用いて下記のよ
うに派生させた．
(2.13)

Pmax.plot [i] ∼ N (Pmax.sp , sdplot )(2.14)

Pmax.prd [i, j] ∼ N (Pmax.plot [i], sdprd )ここで，sdplot および sdprd はそれぞれ試験区固有パラメータおよび調査区間固有パラメータの標準偏差を示し，両パラメータも推定の対象とした．ベイジアンキャリブレーションによるパラメータ推定においては，式（2.3）の Pmax へ試験区 i における調査区間 j に固有のパラメータである Pmax.prd [i, j] を代入して尤度を計算する．モデルを予測に用いる場合には，樹種共通
パラメータである Pmax.sp を代入してシミュレーションを行う．
ベイジアンキャリブレーションによる推定の対象としたパラメータおよびその事前分布について表 2 にまとめた．上記以外のパラメータについては過去の研究に基づいて固定した（表 2）．
ベイジアンキャリブレーションには WinBUGS（Spiegelhalter et al., 1996）を用い，マルコフ連鎖の数を 3 つ，全体の繰り返し回数を 10 万回として MCMC によるサンプリングを行った．初期の 9 万回を除く 1 万回の試行結果から一定区間でサンプルを抜き出し，各パラメータについて事後分布からの擬似標本（以降，MCMC 標本と称す）を得た．

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定
表 3.

315

各パラメータにおける MCMC 標本の概要．抜けている調査区間番号のデータは検証に
用いた．

3. 結果

推定の対象とした各パラメータについて，得られた MCMC 標本における上位，下位四分位点および中央値を表 3 に示す．パラメータ推定が適切に行われたかを判断するため R̂ 統計量（Gelman and Rubin, 1992）によってパラメータの収束判定を行ったところ，ほとんどの変数で判定の基準値である 1 に近いものであったが，大きく超える値となった変数もあった．さらに，各パラメータについて MCMC 標本の時系列変化を確認し（図 2），また MCMC 標本から確率密度を計算して事前分布と比較して（図 3），パラメータ推定における結果の妥当性を検討した．
MCMC 標本の時系列変化に傾向性はみられなかったことから定常状態にあると判断した．事
前分布と比較して MCMC 標本の確率分布が著しく分布の幅が狭くなったり，分布の中心が移動したわけではなかったが，事前分布よりは尖った分布をしていたパラメータが多く，ある程

統計数理

316

図 2.

第 61 巻

第 2 号 2013

各パラメータにおける MCMC 標本の時系列．

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定

図 3.

各パラメータにおける MCMC 標本と事前分布の確率密度の比較．実線は MCMC 標
本，点線は事前分布の確率密度．

317

統計数理

318

図 4.

第 61 巻

第 2 号 2013

Pmax.sp および Pmax.plot に関する MCMC 標本の確率密度．

度は事前分布がデータによって更新されたとみなした．しかし，β5 C と β6 C については事前
分布とほとんど変わらない分布形を示しており，この 2 つの変数については事前分布がデータ
によって更新されていない．本研究においてはベイジアンキャリブレーションによって推定されたパラメータは，事前分布から大きく更新されたものではないものの，妥当性のあるものであると判断して MCMC 標本における中央値を各パラメータの推定値とみなすこととした．幹，枝および根の気温と呼吸速度を表現するモデル式のパラメータ β5 C および β6 C については，下記の検証において出力値を確認したところ呼吸速度は妥当な数値を示したので，他のパラメー
タと同様に取り扱うこととした．
光合成速度の最大値を示すパラメータである Pmax について，スギの樹種平均パラメータと試験区ごとのパラメータの MCMC 標本分布を比較したところ，試験区による違いが顕著であった（図 4）．試験区固有パラメータ標準偏差は 0.143，調査区間固有パラメータの標準偏差は 0.139 と推定され，同じ試験区でも調査区間によって光合成速度を規定するパラメータが大きく変動する可能性が示唆された．
まず，モデルの振る舞いと光合成最大速度パラメータの階層化による効果を確認するため，試験区固有パラメータ，調査区間固有パラメータおよび樹種共通パラメータをそれぞれ用いて，パラメータ推定に用いたデータにモデルをあてはめて林分成長量の再現性を確かめた（図 5（a））．調査区間固有パラメータを用いた場合には観測された林分成長量をよく再現できていた．
試験区固有パラメータを用いた場合にはおおむね良好な推定となったが，顕著な過大推定となるデータが存在した．樹種共通パラメータを用いた場合には，さらに推定誤差が大きくなった．
さらにモデルの信頼性を確認するため，検証用データに樹種共通の光合成最大速度パラメータを用いてモデルをあてはめた（図 5（b））．この検証においてもバイオマス成長量の推定は比較的良好な結果となるデータと，顕著な過大推定となるデータが混在した．
4. 考察

マルコフ連鎖モンテカルロ法によるベイジアンキャリブレーションによって，収穫試験地における継続調査データから林分成長モデルのパラメータ推定を行ったところ，ほとんどの変数
において妥当な推定結果となった（表 3）．しかし，MCMC 標本の確率分布が事前分布から大きく更新されたわけではなかった（図 3）．本研究においては MCMC サンプリングが発散することを避けるため，事前分布の標準偏差を小さく設定した．モデルにおいて光合成による生産量
はポテンシャル生産量に律速係数を掛け合わせる形で計算される（式（2.7））．ポテンシャル生

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定

図 5.

319

モデルシミュレーションによるバイオマス成長量推定の検証結果．

産量を規定する Pmax と温度律速係数パラメータ β1 および β2 がパラメータ推定の対象となっているので，両者が補完的に作用することで MCMC 法によるサンプリングが幅のより狭い分布へと収束しない結果となった可能性が高い．しかし逆に，パラメータ推定が難しいモデル構造に対して，過去の知見を活用できるベイジアンキャリブレーションが有効に機能したと考えることもできるだろう．
試験区固有パラメータは樹種平均的パラメータに対して比較的大きな幅で変動する結果となった（図 4）．7 試験地 13 試験区のデータに対してこれだけ大きな差を生じたことは，遺伝的な影響もあると考えられるが，スギが立地条件の違いに対して明瞭に反応することを示唆している．各試験区における調査区間固有パラメータについても，同じ試験区内においてバラツキが大きい結果となった．本来であれば，試験区の立地的な違いは試験区固有パラメータで代表され，調査区間固有パラメータのバラツキが小さいことが望ましい．期待に反して調査区間固有パラメータのバラツキが大きかった原因の一つとして，加齢によって葉の光合成効率が低下
したり，樹木の水分通導性が悪化したりすることが考えられる．原型となった 3-PG モデルにおいても加齢による律速係数を導入しており（Landsberg and Waring, 1997）樹木生理に関する知見を深めて，加齢の効果をモデルに組み込む必要があるだろう．
林分成長モデルに初期サイズおよび気象値を入力してシミュレーションを行って観測された
林分成長と比較したところ，試験区固有パラメータを用いた推定では多くの調査区間でおおむ
ね良好な結果となったが，顕著な過大推定となる調査区間が存在した（図 5（a））．また，樹種共
通パラメータを用いたシミュレーションによるモデルの検証においても，顕著な過大推定となる場合があった（図 5（b））．過大推定となった調査区間についてデータを確認したところ，ほとんどが 10 年間隔の調査区間もしくは間伐直後の調査区間であった．本研究においては MCMC
サンプリングのシミュレーションにおいて部位別バイオマスを初期サイズに固定して計算を行なっており，シミュレーションの後期になるにつれ実際のバイオマスと初期値との乖離によって誤差が大きくなる．特にシミュレーション期間が長くなる場合には誤差が累積するため，10年間隔の調査区間において誤差が大きくなったと考えられる．間伐後の林分においては，間伐によって林冠に孔隙が生じて葉群を通過しないことで林分としての日射吸収様式が変化する．
本研究においては樹種で固定の吸光係数を用いるため，このような変化に対応できず，実際よ

統計数理

320

第 61 巻

第 2 号 2013

りも多く日射量を吸収したように計算されることによって過大推定となったと考えられる．このような状況に対応するために計算方法の改善や，間伐の効果などを反映できるようにモデル
構造の改良などに取り組む必要がある．
本研究においては，収穫試験地という林業の生産量予測を目的としたデータを用いて，森林の成長過程を気象環境への応答として記述する林分成長モデルのパラメータ推定を行うことができた．その過程においては，パターン試行モデリングの考え方（山道・角谷, 2009）およびその実行手段としてのベイジアンキャリブレーションが有効であった．本研究で開発したモデルは林分構造を設定すれば，全天日射量，平均気温，最高気温および相対湿度の月別値という入手が容易な気象値さえあれば適用することが可能であり，平均的なスギ人工林の動態が予測できる．擬似的な気象値を与えれば，気候変動がスギ人工林における木材供給サービスにどのような影響を与えるのかをシミュレーションすることもできる．さらに，限定的ではあるが立地条件等によるバラツキを示すパラメータ（sdplot）と樹種平均的なパラメータ（Pmax.sp）を利用
してランダムにパラメータを生成する繰り返しシミュレーションが可能となり，不確実性を考慮した予測が可能となる．このように本モデルは環境リスクの評価のための可用性が高いツールであると言えるだろう．収穫試験地データのような時系列データは貴重なものであるが，古くからの計測であるが故に計測項目が少ないといった弱点もある．しかし，そのようなデータであっても，パターン試行モデリングの考え方に基づいて環境リスクを評価できるようなモデルを設計し，ベイジアンキャリブレーションを応用することで上手くパラメータ推定を行うことができることが示された．このように収穫試験地の継続調査データはより幅広く利用できるようになると期待でき，有効活用を推進するためにもさらに調査を継続していくことが肝要で
ある．
謝

辞

これまで収穫試験地の管理，調査に関わった全ての国有林および森林総合研究所職員のみなさまに対して，ここに深甚なる敬意と謝意を表します．また，本研究は農林水産技術会議委託
プロジェクト研究「農林水産分野における地球温暖化対策のための緩和及び適応技術の開発」によって行われた．

参 考 文 献Aber, J. D. and Federer, C. A.（1992）. A generalized, lumped parameter model of photosynthesis,
evaporation and net primary production in temperate and boreal forest ecosystems, Oecologia,
92, 463–474.
Bergh, J., Freeman, M., Sigurdsson, B., Kellomäki, S., Laitinen, K., Niinistö, S., Peltola, H. and Linder,
S.（2003）. Modelling the short-term eﬀects of climate change on the productivity of selected
tree species in Nordic countries, Forest Ecology and Management, 183, 327–340.

千葉幸弘（2004）. CO2 収支評価における森林構造と動態，地球環境，9，213–220.
Coops, N. C., Waring, R. H. and Law, B. E.（2005）. Assessing the past and future distribution andproductivity of ponderosa pine in the paciﬁc northwest using a process model, 3-PG, EcologicalModelling, 183, 107–124.
Fukuda, M., Iehara, T. and Matsumoto, M.（2003）. Carbon stock estimates for sugi and hinoki forests
in Japan. Forest Ecology and Management, 184, 1–16.
Gelman, A. and Rubin, D. B.（1992）. Inference from iterative simulation using multiple sequences,
Statistical Science, 7, 457–511.

人工林長期継続調査データを利用した林分成長モデルのパラメータ推定321

細田和男，光田 靖，家原敏郎（2010）. 現行立木幹材積表と材積式による計算値との相違およびその修正方法，森林計画学会誌，44，23–40.
細田和男，西園朋広，佐野真琴，高橋與明，光田 靖，家原敏郎（2012）. 間伐強度の異なるヒノキ人
工林における 35 年間の炭素固定量の比較，森林計画学会誌，46, 1–14.
Landsberg, J.（2003）. Modelling forest ecosystems: State of the art, challenges and future directions,
Canadian Journal of Forest Research, 33, 385–397.
Landsberg, J. J. and Waring, R. H.（1997）. A generalised model of forest productivity using simpliﬁed concepts of radiation-use eﬃciency, carbon balance and partitioning, Forest Ecology andManagement, 95, 209–228.
松本光朗（2008）. 地球温暖化と森林，森林科学, 52, 4–8.
松本陽介，重永英年，三浦 覚，長倉淳子，垰田 宏（2006）. 温暖化に対するスギ人工林の脆弱性マップ，地球環境, 11, 43–48.
Mitsuda, Y., Hosoda, K., Iehara, T. and Matsumoto, M.（2010）. Parameterization of a process-based
forest growth model using long-term yield survey plot data for predicting carbon sequestrationin Cryptomeria japonica plantations, Carbon Sequestration: Methods, Modeling and Impacts（eds. E. Hoch and S. Grunwald）, 175–202, Nova Science Publishers Inc, New York, USA.
Monsi, M. and Saeki, T.（1953）. Über den Lichtfaktor in den Pﬂanzengesellschaften und seine Bedeutung für die Stoﬀproduktion, Japanese Journal of Botany, 14, 22–52.
Mori, S., Yamaji, K., Ishida, A., Prokushkin, S. G., Masyagina, O. V., Hagihara, A., Hoque, A. T. M. R.,
Suwa, R., Osawa, A., Nishizono, T., Ueda, T., Kinjo, M., Miyagi, T., Kajimoto, T., Koike,T., Matsuura, Y., Toma, T., Zyryanova, O. A., Abaimov, A. P., Awaya, Y., Araki, M. G.,
Kawasaki, T., Chiba, Y. and Umari, M.（2010）. Mixed-power scaling of whole-plant respiration from seedlings to giant trees, Proceedings of the National Academy of Sciences, 107,
1447 –1451.
農林水産技術会議事務局（2011）. 地球温暖化が農林水産業に及ぼす影響評価と緩和及び適応技術の開発，農林水産技術会議研究成果，No. 483，農林水産技術会議事務局．
大谷義一（2008）. 生態系炭素の観測体制，森林科学，52，19–22.
Pretzsch, H., Grote, R., Reineking, B., Rötzer, Th. and Seifert, St.（2008）. Models for Forest Ecosystem Management: A European Perspective, Annals of Botany, 101, 1065 –1087.
志賀弘行（2003）. アメダス観測地点における最低気温を用いた水蒸気圧の推定，北海道立農業試験場
集報，84，99–100.
森林総合研究所（1996）. 森林長期モニタリングシステム 収穫試験地の時系列データの収集と整備 ，
森林総合研究所，つくば．
Spiegelhalter, D., Thomas, A., Best, N. and Gilks, W.（1996）. BUGS: Bayesian Inference Using
Gibbs Sampling, Version 0.50, MRC Biostatistics Unit, Cambridge.
山道真人，角谷 拓（2009）. マルコフ連鎖モンテカルロ（MCMC）法を用いたシミュレーションモデル
のパラメータ推定：ベイジアンキャリブレーション入門，日本生態学会誌，59，207–216．

322

Proceedings of the Institute of Statistical Mathematics Vol. 61, No. 2, 307–322(2013)

Parameterization of a Forest Stand Growth Model Using Long-term Field
Survey Plot Data
Yasushi Mitsuda1 , Kazuo Hosoda2 and Toshiro Iehara2
1 Shikoku Research Center, Forestry and Forest Products Research Institute;

Now at Faculty of Agriculture, University of Miyazaki
2 Department of Forest Management, Forestry and Forest Products Research InstituteWe developed a matter balance based stand-level forest growth model for Cryptomeria japonica planted forests using long-term field survey plot data. This model consists ofsix processes: (1) photosynthetically active radiation absorption; (2) conversion to grossprimary production; (3) constraints on photosynthesis by environmental factors; (4) respiration; (5) litterfall and root turnover; and (6) biomass partitioning. For parameterization,
we estimated time-series biomass and biomass growth data using repeated measurementsof long-term field survey plots. We applied Bayesian calibration methodology to parameterization in our process-based model using time-series biomass and biomass growth data
and time-series climatic values of solar radiation, mean temperature, and vapor pressure
deficit. The results of model validation indicate that, although some large over-estimation
errors were observed, our model could represent patterns of biomass growth measured inlong-term field survey plots.

Key words: Pattern-oriented modeling, Bayesian calibration, Cryptomeria japonica planted forest.