© オペレーションズ・リサーチ

森林資源管理における最適化モデルの展開 (II)
―森林単位での最適化モデル―
吉本

敦，木島

真志

前稿では，経営体の最小単位となる 林分単位での最適化モデル について詳説した．本稿では，複数の林分からなる森林単位の最適化モデルについて，これまで展開されてきた古典から近年まで，その関連研究について
著者の視点から詳説する．

キーワード：伐採計画最適化，線形計画法，混合 0-1 整数計画法

1. 古典的な森林単位の問題：許容伐採量の探索がある．
面積・材積を制御する折衷法は上記二つを組み合わ

歴史的に多くの注目を集めていたもう一つの問題は，せたすべてのものを含む．それゆえ，必要とする情報

森林単位において許容伐採量 (Allowable Cut: AC)

如何によりさまざまな方法が考案できる．この方法のを如何に決定するかというものである．鈴木 [1] の収特徴は，最終的な解がどちらの基準も満たすということ

穫予定論にさまざまな決定法が示されている．大別すである．たとえば，得られた解が現存の森林の伐採に必ると三つに分かれる．一つが森林面積の制御による区要な計画期間が希望する期間と一致するかどうか，そし画輪伐法 (Area Control)，二つ目が材積配分法などのて各期の伐採面積が許容できるかどうかなどがチェッ

材積の制御による決定法 (Volume Control)，そして，

クされ，もしそれらの条件を満たさなければ，再度解その折衷方法 (Volume-Area Check) である．

の探求が行われる．したがって，コンピュータプログ

区画輪伐法は，与えられた森林を 1 回の伐期内で，各期同じ面積を伐採することにより，現存の森林を制御されたものにすることを目的としている．均一の林

ラミングによる解の探求が必要不可欠となる．

2. 確率過程の応用と収穫予測分からなる森林であれば，許容伐採量の決定と伐採の古典的なアプローチとは別に国レベルの収穫予測を

配分は簡単に行うことができる．伐期齢を決定すれば行うために確率過程を用いたアプローチが 1960 年代

各期で伐採される面積は総面積を伐期齢で割れば簡単

初期に発表された．当時，科学技術庁資源局よりわがに求められる．この方法では林分の伐採ルールとして，国の木材収穫を予測するマクロ的な課題が林学研究者「最も古い林分から伐採する」が採用される．したがっ

に問われた．そこで鈴木 [2, 3] が世界で初めて確率過て，それぞれの林分の林齢と材積がわかれば，どの林程を駆使した方法を提示した．鈴木 [4] 氏は，齢級ご分をいつ伐採し，各期どれだけの伐採量があるかを計とで伐採される割合を 減反率 と定義し，マルコフ

算することができる．

連鎖の出生死滅過程 (birth–death process) を応用す

材積を制御する方法では，計画期間での伐採量を一ることにより，各齢級における伐採−成長の過程を表定にするように計画が立てられる．この種の方法には現し，将来的な林齢の遷移と収穫の予測が可能になる

さまざまな公式が開発された．その代表的なものには．なお，齢級とは林齢を単位年数ことを示した（図 1）

材積配分法をはじめ，Grosenbaugh 式，Hundeshagen（日本では 5 年，アメリカでは 10 年）でくくった森林

式，Von Mantel 法，Hanzlik 式，オーストリア式などの年齢を表現したものである．

よしもと あつし
統計数理研究所〒 190–8562 東京都立川市緑町 10–3

yoshimoa@ism.ac.jp
このしま まさし
琉球大学農学部〒 903–0213 沖縄県中頭郡西原町字千原 1 番地

konoshim@agr.u-ryukyu.ac.jp

図1

単純マルコフ連鎖による伐採プロセス

608（42）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ

第 t 期において，第 i 林齢の森林は次の期には確率

pti,i+1 で第 (i + 1) 林齢に遷移する．すなわち，

ば，単純マルコフ連鎖の下で，林齢の遷移確率行列を用いて，将来的な木材供給量の予測が可能になる訳である．
ここで問題となるのは，各林齢の遷移確率行列をど

pti,i+1 = Pr{yt+1 = i + 1 | yt = i}のように推定するかである．今，森林が植林され，丁度yt は第 t 期の林齢状態を表す変数である．したがって，

第 i 齢級で伐採される，すわなち第 i 林齢から第 (i + 1)

その補完の確率 (1 − pti,i+1 ) は，第 t 期から第 (t + 1)

林齢になるまでに伐採される確率を第 i 齢級の減反率

期に遷移する間に第 i 林齢の森林が伐採される確率とqi と定義する．その結果，第 i 林齢の森林 (yt = i) がなる．なお，ここでは期間間隔と林齢間隔，また齢級第 (i + 1) 林齢になる前に伐採される確率は下記の条の林齢の単位は同じとしている．また，伐採後は即座件付き確率で表され，この確率が上記林齢の遷移確率に再造林されるものとし，その状態を第 1 林齢と表す．

行列の要素となる．

このように yt の遷移過程を定義すると図 1 のような

pti,1 = Pr{yt+1 = 1 | yt = i}

単純マルコフ連鎖により林齢の遷移の様子を表現することができる．

=

次に (n × n) の林齢の遷移確率行列を定義する．こ=の行列は林齢状態の遷移を描写するものであり，それ

⎜
⎜ pt
⎜ 2,1
⎜
t
Pt = ⎜⎜ p3,1
⎜
⎜ ·
⎝
1

qk
k=1

かを示すものである．

pt1,1

(4)

i−1

1−

ぞれの状態から次の状態にどれだけの割合が遷移する⎛

Pr{yt+1 = 1, yt = i}Pr{yt = i}
qi

pt1,2

0

0

···

0

pt2,3

0

···

0

0

pt3,4

···

·

·

·

···

0

0

0

···

⎞

ここで減反率と区別するために，この確率を第 i 齢級

⎟
0 ⎟
⎟
⎟
0 ⎟
⎟ (1)
⎟
· ⎟
⎠
0

の伐採率としている．減反率が植林後それぞれの齢級

0で森林が伐採される確率で，伐採率は植林後の各林齢まで成長し，次の林齢になる前に伐採される確率である．鈴木氏は樹木の直径成長を単純なポアソン過程で表し，次に森林所有者の伐採行動を 直径がある望む径に達したら伐採する とし，この待ち時間が，植林

この林齢の遷移確率行列を用いると，現状がわかれば，されてから伐採されるまでの時間となり，ある齢級に将来的な林齢の状態が下記により予測可能となる．

おいて望む径に成長し伐採される割合，すなわち 減

Xt+1 = Xt · Pt t+1 t+1 t+1

x1 x2 x3 · · · xt+1
= xt1 xt2 xt3 · · · xtn
n
⎞
⎛
0
··· 0
pt1,1 pt1,2 0⎟⎜⎜ pt
pt2,3 0··· 0 ⎟⎟
⎜ 2,1 0
⎟
⎜
t
t
⎟ (2)
·⎜
p
0
0
p
·
·
·
0
3,4
⎟
⎜ 3,1
⎟
⎜
⎟
⎜ ·
·
·
·
·
·
·
·
⎠
⎝
1
0
0
0
··· 0

反率

がガンマ分布から算出できるとした．

実際，わが国では，全国森林計画に即して各地域で作成される地域森林計画において，ガンマ分布を仮定
した減反率分布と地域の齢級別森林面積を示した森林
資源構成表から，長期の森林資源量の変動を予測している [5]．藤掛 [6] はガンマ分布の推定法の最尤法による改良を試みた．しかし，根本的に実際の伐採傾向とガンマ分布に基づく減反率からなる伐採傾向に大きな乖離があることが指摘された [7, 8]．その後，減反率

ここで Xt は第 t 期の各林齢ごとの森林面積ベクトルを推定するモデルの改良が Yoshimoto [9–11] によりを表し， は転置を示す．仮に各林齢に対する単位面積行われた．いずれにせよ，今日ではポアソン過程−待当たりの収穫材積がわかれば，木材供給量 Vt の予測もち時間問題，マルコフ連鎖による出生死滅過程は一般

以下の式により可能になる．

的な確率過程の教科書では当たり前のように出てくる

n

Vt =

例題ではあるが，当時，鈴木氏がそれらを組合せ林学

vi · xti · pti,1(3)への応用の道を開拓したという点はほかに例を見ない．

ただし，vt は第 i 林齢における単位面積当たりの材積を

3. 森林単位での最適化モデル（線形計画法・混合 0-1 整数計画法の応用）

i=1

示す．このように，現在の各林齢ごとの森林面積とその伐採される確率，さらに単位面積当たりの材積がわかれ

2022 年 11 月号

さて，70 年代に入ると，許容伐採量の探求ととも

Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited. （43）
609

表1

計画期間 10 期最小伐期齢 6 期の実行可能施業案に伐採計画の最適化問題について線形計画法の適用が画は 5 年ごと，あるいは 10 年ごとを単位とするため，本格的に展開され始めた．伐採計画問題に対する線形

毎年の単位ではなく，1 期＝ 5 年または 10 年として計画法の最初の応用については，諸説あると思われる

計画を立てる．このような意味合いからも，線形計画

が [12]，なかでも Johnson and Scheurman [13] が提法による計画は長期的な方向性を示すものとして捉え案した Model I と Model II は今日の伐採計画問題のられていた．表 1 で 1 回だけ伐採される施業案が 1，基本となっている [14]．当初はコンピュータのパワー

6，10，13，15〜20 であり，それぞれ第 1 期目から順不足で対象とする数百あるいは数千にも及ぶ林分から

次 10 期目の伐採に対応する．最小伐期齢を 6 期とすなる森林に対し，森林を樹種，樹齢，林道などの類似状

れば，計画期間内では最大 2 回の伐採が可能となる．

態によって区分層 (Strata) に分類し，問題規模の縮小

上記以外の施業番号が 2 回の伐採を伴う施業案であ化を図っていた．Model I と Model II の違いは同じ区

る．それぞれの区分層に対し，考えられる施業案を構

分層あるいは林分に対して時期の異なる複数の伐採を築し，一つ一つの施業案に対応する森林面積が決定変どのように捉えるかにある．計画期間内で 1 回しか伐数となる．

採を考慮しなければこれらのモデルに違いはない．定決定変数 {xij ：第 i 番目の区分層における第 j 番目の式化に当たっては，グルーピングされたそれぞれの区

施業案適用の森林面積割合}

分層に対し，いつどれだけの面積を伐採するのかを決

決定変数が決まれば，次は目的関数において決定変定するため，それに対応するような決定変数を考える

数に対応する係数 {cij } を推定する．たとえば収益の必要がある．そこで，Model I では与えられた計画期最大化であれば，それぞれの施業案から計画期間内で間内で考え得る伐採パターンを想定し，それぞれのパ得られるであろう単位面積当たりの収益を算出すればターンに対して，どれだけの面積を割り当てるかを決よい．単純に考えれば，まず各区分層における材積成

めることにより，全体の問題を定式化する．それぞれ

長曲線から伐採が可能な期における ha 当たりの材積のパターンを各区分層に対する施業案 (treatment) とを算出する．次に立木から丸太への歩留まりを掛ける

すると，たとえば，計画期間 10 期の場合，最小伐期齢

ことにより，搬出可能な材積を推定し，立米当たりの丸（伐採を行える最小の齢級）を 6 期とすれば，表 1 のよ

太価格を掛けることにより，各期での伐採から得られ

うに 20 の伐採のパターンが考えられる．表内の 0 はる ha 当たりの収益がだいたい算出できる．したがっ何もしないことを意味し， 1 は現存の森林の伐採，そて，複数回の伐採に対しては各伐採に対応する収益をして 1 より大きな数字は再造林後の伐採を意味し，そ

合計することにより，決定変数の係数が算出でき，目の番号は伐採される齢級を示す．なお，一般に伐採計的関数が下記のように特定できる．ただし，N は区分

610（44）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ

図2

Model II における伐採連鎖のネットワーク (Yoshimoto and Asante [15])層の数を示す．区分層の現時点における樹齢が最小伐の制約がある．これは，少なくとも計画期後の将来的期齢に達していない場合は，係数 {cij } を 0 にすればな経営の再開を可能とするような制約と考えることがよい．

できる．また，最後の計画期における林齢分布に対すNる制約も考慮できる．現在は昭和 30 年代の拡大造林

20

cij · xij

max

(5)

の影響もあり，日本全体の人工林の林齢分布が 50 年生

目的関数が定義されれば，次は制約条件を考える必

多様性の観点から見れば，あまり好ましい状態とはい要がある．一般的な伐採計画では計画期間内での森林

えない．それよりは若齢林から老齢林まで満遍なく森の売買は想定しないため，森林面積に対する制約が発林が分布している方が多様性があるともいえる．そう

生する．すなわち，現在の与えられた森林に対して計した条件に対してもここで定義した決定変数を用いて画を組む．決定変数が各施業実施の森林面積であるの定式化が可能である．なお，Model I では，一つのマで，各区分層に対応する決定変数の総和がその区分層

ネジメントユニットである区分層に対し計画期間内での総面積以下であればよい．

考えられる施業案を想定し，それぞれの施業案に対し

i=1 j=1

あるいは 60 年生に偏っている．こうした偏りは，生物

目的を達成するように伐採計画が行われる．したがっ20

xij ≤ Li

∀i

(6)

j=1

て，計画期間内では区分層内の森林が他の区分層に移行することは考えられていない．

ただし，Li は第 i 番目の区分層の森林面積である．

一方，Model II では，計画期間内で一部の森林の次に採用される制約には，計画期間を通した伐採量区分層から他の区分層に移行することも考慮できる．

に関するものがある．雇用の持続性などを考慮すれば，

Model II の決定変数は各期における伐採施業面積とし，期間ごとの伐採量の急激な増減は望ましくないという同一区分層における複数の伐採については，伐採連鎖の理由から伐採量はなるべく一定あるいは減少しないとネットワークを使用し規制している．図 2 に Model II

いう制約が課せられる．各施業案において， 1 回目ので考慮される伐採連鎖のネットワークを示す．

伐採かそれ以降の伐採かで伐採量は異なるが，材積のノード 0 が始点でノード（T + 1）が終点となる．

成長予測により，伐採量の算出が可能になる．その結その他の番号は伐期を表す．伐採の連鎖を表現するたt
を t 期における i 番目の区分層で j
果，たとえば vij

めに，一連の植林−伐採を考え，第 i 番目の区分層番目の施業を行って得られる伐採量とすれば，i 番目

で s 期に植林された森林が t 期に伐採される決定変の区分層から t 期に得られる伐採量 Vit は下記で計算数を xist とする．ただし，伐採は t = T まで有効で，できる．

t = T + 1 の場合，T 期まで伐採されないことを意
Vit =

20

味する．たとえば i 番目の区分層に対し現状の森林を
t
vij
· xij

(7)

j=1

後 1 期目で植林され T 期目で伐採する施業は xi1T とその結果，伐採量に関する制約は下記のようになる．
N

1 期目で伐採する施業の決定変数は xi01 となり，そのなる．このような伐採の連鎖を規制することにより複数の伐採を取り扱うことができる．xi0t は現状の森林

Vit = v0 ,t = 1, . . . , T

(8)

i=1

が t 期で伐採される施業に対する決定変数である．したがって，ある区分層に対し，下記の制約により一連

ただし T は計画期間で，v0 は一定伐採量である．そのの伐採の連鎖を制御できる．なお，k は最小伐期齢で他の制約としては計画期の最後の状態に対する資源量ある．

2022 年 11 月号

Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited. （45）
611

表2

Model I と Model II の定式

注 N : 区分層数，H: 施業数，T : 計画期間，k: 最小伐期齢

max(0,s−k)xiτ s =

τ =0

れる伐採計画最適化モデルを開発し [16]，さまざまな

T +1xist , ∀i, s = 1, . . . , T

t=min(s+k,T +1)

(9)
T +1

xi0t = Li , ∀i

(10)

t=1

デルで，水源地の保全や野生動物の保護面積などを制約条件として導入することで，木材生産といった経済活動が環境保護に及ぼす影響について線形計画法の枠伐採量に関する制約についても下記で対応できる．
N max(0,t−k)i=1

社会ニーズに対応するようになった． FORPLAN は
基本的にソルバーのインプットファイルを生成するモ組みで定量的に評価することを可能とした．
Ni
vst· xist = v0 ,t = 1, . . . , T

(11)

N

xi0 ≥ L̃ :Model I :Model II :

i=1

s=0

(13)

i
ここで vstは i 番目の区分層が s 期で植林され t 期で伐採されるときの伐採量である．Model I 同様に他のxi0(T +1) ≥ L̃i=1

L̃ は保護面積，xi0 ，xi0(T +1) はそれぞれ Model I と制約もここで定義する決定変数により対応できる．目

Model II での無伐採面積である．上記のように環境問的関数は下記のように定式化する．

題を木材生産に対する制約として捉えるのではなく，そN T +1 max(0,t−k)max
i=1 t=1

れら森林の保有する環境的な機能も同時に最適化の一部

cist · xist

(12)

s=0

とする多目的計画法の応用研究も展開されている [17]．
区分層を利用した FORPLAN が実際の森林管理に

ここで cist は i 番目の区分層が s 期で植林され t 期で利用されるようになってから，伐採を担う現場から，伐採される際の貢献あるいは収益の現在価値などであ「どの場所の林分を伐採するのか」といった疑問が浮上

る．表 2 に Model I と Model II の定式を示す．

してきた．そこで線形計画法から混合 0-1 整数計画法

森林資源を取り巻く社会ニーズは，これまでの伐採への移行が始まった．Model I，Model II の伐採施業

量の維持による地域雇用の確保だけでなく，環境保全に関わる決定変数，係数を i 番目の区分層から i 番目を訴える傾向が出てきた．米国では 1969 年に制定さの林分としバイナリー変数に置き換え，伐採量制約をれた国家環境政策法 (National Environmental Policy

変動許容率 α を用いて調整すれば表 3 のように離散型Act) に基づき国有林の森林計画策定段階において環境

Model I と Model II が構築できる．

アセスメントが求められるようになった．そこで，アメリカの国有林を管理する森林局 (Forest Service) は

3.1 混合 0-1 整数計画法による隣接問題への対応(Unit Restriction Model)

管理計画の策定に際して，上記 Model I と Model II

80 年代後半に入りマダラフクロウ保護などさまざまを採用した FORPLAN (FORest PLANning) と呼ばな野生動物保護問題が浮上した．その中で注目された

612（46）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ

表3

離散型 Model I と Model II の定式

注 N : 林分数，H: 施業数，T : 計画期間，k: 最小伐期齢のが大面積皆伐や森林に生息する野生動物の生息地のきない．(14) 式はその制約に対応している．Model II

分断・劣化の回避を狙った隣接制約 (adjacency con-においては，緑化再生期間を g とし，下記となる．

straint) であった．後に，この種の制約を扱ったモデル

max(0,t−k)を Unit Restriction Model (URM) と呼ぶようになった [18]．この制約では， 隣接し合う林分同士の同時期
伐採は不可 とした．基本的に隣接するもの同士の伐

s=0xist +(t+g−1) max(0,t+g−1−k)p=t+1xjsp ≤ 1,

s=0∀j ∈ NBi , ∀i, t = 1, . . . , T − g + 1

採を回避する隣接対 (pairwise) の制約式 (xi +xj ≤ 1)

(15)で十分であったが，当時は PC の容量も少なく，それ3.2 混合 0-1 整数計画法による集約化問題への対ゆえ 如何に少なく制約式を設定できるか の議論が

応 (Area Restriction Model)．この制約の延長が緑化再生進んだ（文献 [19] など）

隣接制約の導入が始まった後，実際の現場で施業費

制約 (green-up constraint) というもので隣接の林分

用が増す小規模な林分への過度な分断が発生し，費用の伐採の遅延により伐採箇所の緑化再生を促すもので削減への規制緩和が始まった．そこで，隣接制約に対ある．伐採の遅延の期間が 1 期間であれば，隣接制約

し，林分の集約をある程度許容する最大許容伐採面積規と同じになる．表 4 に隣接制約を考慮したそれぞれの制 (maximum opening size requirements) が導入さモデルを示す．Model I では施業隣接行列 (activity

れた．この規制はある程度の規模面積までの隣接伐採

adjacency matrix) [20] という施業間で緑化再生期間を許容するもので，隣接制約が分断化 (fragmentation)の同時期伐採を (0,1) の要素 avhl で示した行列を導入を促すのに対し，これは集約化 (aggregation) も許容

し，下記制約で対応する．

するものである．日本においても低コスト林業の推進

avhl · xih + avlh · xjl ≤ 1, ∀j ∈ NBi , ∀i, ∀l(= h), ∀h(14)に伴い施業費用の削減に向けた林分の 団地化 が進められている．URM に対し，この種のモデルは Area

Restriction Model (ARM) と呼ばれている．

なお，NBi は i 番目の林分に隣接する林分のインデック

集約化に伴う制約式の構築については，ルールに基スセットである．ここで avhl が 1 の場合は h 番目の施づき候補となるクラスターの生成をもとにしたものが業と l 番目の施業が同時期の伐採を行うことを意味し主流であった．PATH 定式 [21]，Cluster 定式 [22]，ている．その場合，i 番目の林分に h の施業が施される

Bucket 定式 [23] がそれらである．集約化は境界線あると，隣接する j 番目の林分では l の施業を施すことがで

いは，境界点を共有する林分同士に加えて，それらに対2022 年 11 月号

Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited. （47）
613

表4

隣接制約と離散型 Model I と Model II の定式

注 N : 林分数，H: 施業数，T : 計画期間，k: 最小伐期齢

図3

森林マップ例，a) 空間配置，b) グラフによる描写

してさらに逐次的に連続した林分の集合体を生成してる必要があるため，対象とする林分の数，計画期間のいく必要がある．吉本ら [24] や Yoshimoto et al [25]

サイズに制限があった．それに対し，Yoshimoto andは，予め望む集約形状を考慮し，まず，逐次的な連結操Asante [26] は，最大フロー問題を用いて集約化を捉え作のルール設定に基づき林分集合体の候補を構築する

るアプローチを提示した．

方法を提示した．そして，それら候補の中から，重複

図 3 に森林マップの例を示す．ここで円内の番号はする林分集合体同士の同時期選択を回避することによ

林分の ID，その下の数値が面積を示す．a) の空間配り，最適な林分集合体の組み合わせを探索するとした．

置をノードとアークからなるグラフで表すと b) とな

なお，重複の回避は隣接制約に類似するものである．

る．ここではスーパーノードとして架空のノード 0 を

上記の方法は事前に候補となるクラスターを準備す

追加し，すべてのノードからの連結を仮定する．グレー614（48）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ

図4

最大フローによる集約 a) 三つの隣接，b) 四つの隣接 (Yoshimoto and Asante [26])のポリゴンは保護対象林分であるが，ここでは通常のに，開発されたのが逐次三角連結 (Sequential Trian-

林分として取り扱う．それぞれの面積を各ノードへの

gle Connection: STC) である．これは三つの隣接しフロー，そしてすべてが他のノードを経由するか直接合う林分の結合において，二つの結合が最大フローでスーパーノード 0 に流れ，回収されるものとする．他

選択されたら，選択されていない三つ目のアークを順のノードを経由する場合は，ノード同士の集約でクラ次結合するものである．これまで隣接関係は境界線を

スターを形成し，ほかからの流入なしに直接スーパー基準にしていたが，STC では境界点での隣接関係もノードに流れるものは単独林分となる．次に各ノード利用する．STC のアーク連結には新たにバイナリー

からスーパーノードに流入する量が最大許容面積以下

変数 yijになるようにすれば，単なる集約だけの問題は解決でプロセスで (i, j) 間のアークが連結されれば 1 とす

きる．次に集約に関わる制約を示す．

る．なお，yij = yij

B(n)を導入し，第 (n − 1) 番目の STC 連結
B(1)である．たとえば図 4a) で⎞は二つの隣接 (yik = yjk = 1) が最大フローで結合⎟⎜⎟⎜
yij = 1, ∀i
yi0 +
⎟
⎜
⎟
⎜
j∈NBi
⎟
⎜
⎟
⎜
⎟
⎜
w
ij ≤ L̂ · yij , ∀j ∈ NBi , ∀i⎟⎜⎟⎜
⎜ wi0 +w
w
ij =ji + Li , ∀i ⎟
⎟
⎜
⎟
⎜
j∈NBi
j∈NBi
⎟
⎜
m
⎟
⎜
⎟
⎜
w
=
L̂
i0
⎟
⎜
⎟
⎜
i=1
⎟
⎜
⎟
⎜
wi0 ≤ Lmax⎠⎝yij ∈ {0, 1}, ∀j ∈ N Bi , ∀i

されたとする．その場合，STC の初期条件として，⎛

yij + yji ≤ 1, ∀j ∈ NBi , ∀i

B(1)yik = yjk = yik
B(2)

合 yij

= 1 とし，3 番目の結= 1 が STC により決定される．また図 4b)では (i, l, k) が境界線と境界点で互いに隣接関係にあ

(16)

B(1)

るため，yil = ykl = yil
を導く．同様に yjk =
B(3)

yij

B(0)

2 · uijk + 3 ·積である．ここで注意する点は，最大フローで得られるクラスター内の林分同士は最小のアーク連結で生成されていることである．つまりクラスター内の林分においてすべての隣接するアークが結合している訳ではない．仮にこの制約式群を伐採計画に適用しようとすれば，他の隣接制約により，たとえば三つの隣接し合う
林分からなるクラスターは生成されず，隣同士のみの連結による細長いクラスターのみが生成されてしまう．
最大フローにより得られるクラスター内の林分群
ですべての隣接を結合し，隣接制約を回避するため=1= 1，さらに= yij + yji , ∀j ∈ NBi , ∀i, (i < j)
n−1

(n)

n

(l)

l=1

B(l)

yij

uijk ≤(n)≤ 2 · uijk + 3 ·の結合（林分の集約）を上記の制約で規制でき，最終的にスーパーノードに流れる流量を最大許容面積 Lmax

=

B(2)
yik

を式で表したものが下記である．

の結合を表し，1 なら結合となる． wij は非負の変数で以下に制約すればよいことがわかる．ここで L̂ は総面

B(1)
yjk

B(2)

= 1 が yik

すべての隣接が結合される．このような STC の操作

yiji から j ノードへの流量を示す．流れる流量とアーク

B(1)

= ykl

= 1 となり，本クラスター内の林分群においてバイナリー変数 yij は i 番と j 番目の林分間のアーク

2022 年 11 月号

B(1)

= yjk

B(l)

+ yik

B(l)

+ yjk

l=1
n−1
(l)
uijk + 1,
l=1∀j ∈ MBi , ∀i, ∀k ∈ (MBi ∩ MBj ),n = 1, . . . , N − 1
(n)

n+1

B(l)

yij

3 · uijk ≤

B(l)

+ yik

B(l)

+ yjk

,

(17)

l=1

∀j ∈ MBi , ∀i, ∀k ∈ (MBi ∩ MBj ),n = 1, . . . , N − 1
B(n+1)yij≤(n)uijk ,k∈(MBi ∩MBj )∀j ∈ MBi , ∀i, n = 1, . . . , N − 1
N

B(l)

yij

≤ 1, ∀j ∈ MBi , ∀il=1

Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited. （49）
615

(n)

ここで uijk は第 n 回目の STC プロセスにより互いにものを MF-Model I，Yoshimoto and Asante [15] で隣接する (i, j, k) からなる三角連結の成立を示すバイ

Model II によるものを MF-Model II と命名している．

ナリー変数であり，MBi は i 番目の林分に境界線およまた，Model I に対して施業の概念を使用せずに同時び境界点を通して隣接となる林分のインデックスセッ

期のクラスター伐採を実現することも可能である [15]．
上記の伐採計画に伴う集約と同時に 1990 年代からトである．
以上，集約によるクラスターの形成は最大フローお2000 年代に入ると森林の空間構造の重要性が議論されよび STC により可能になる．最後に空間的な集約と始めた [27, 28]．たとえば，前述のマダラフクロウは，伐採活動を紐づける必要がある．すなわち，Model I

老齢林のような複雑な森林構造をもつ植生が，ある一ではクラスター内の林分に対して同一施業を施すこと

定の面積で空間的に広がった森林に生息するといわれにより，クラスター内すべての林分を同時期に伐採すている [29]．また，野生動植物の生息地を繋ぎ，それるらの移動を可能にするコリドー（緑の回廊）の設計もクラスター伐採 を実現させる．それに対して，

Model II では同時期の クラスター伐採 を実現する

連続性を考慮しなければならない．このような問題をように制約すればよい．Model I を用いると下記の制

Natural reserve network 問題とし，数多くの論文が約が必要となる．

OR 関連の論文誌で発表されているが，基本的にここ

⎛

⎞

h
h
2 · zij≤ xih + xjh ≤ 2 · zij
+ 1,
⎟
⎜
⎟
⎜
∀j ∈ NBi , ∀i, ∀h(= 0)
⎟
⎜
⎟
⎜ N
H
⎟
⎜
B(n)
h
⎟
⎜
yij
=
zij ,
⎟
⎜
⎟
⎜ n=1
h=1
⎟
⎜
⎟
⎜
∀j ∈ NBi , ∀i
⎟
⎜
⎟
⎜
v
v
⎟
⎜ a · x + a · x ≤ 1,jl
lh⎠
⎝ hl ih∀j ∈ NBi , ∀i, ∀l(= h), ∀hで提示した集約問題で同様に対応できる．そのほかにも，上記アプローチはある保護区の拡張を狙った周辺
バッファーとの結合にも対応できる [30]．

(18)

4. おわりに
森林資源の成長において， 時間は非常に重要な投入要素であり，森林は長期的な視野に基づいて管理する必要がある．それゆえ，紀元前から収穫といったh
ここで新たなバイナリー変数 zijを導入し，隣接する

森林資源利用に対する管理が重要視されてきた．そしi と j 番目の林分に h 番目の施業が実施されるときにて，16–17 世紀には，木材供給の予測や輪伐（ローテー

1 となる．施業が異なるが，同一時期での伐採を伴うション）といった，長期的な視野に基づく持続的な資施業に対しては，3 番目の制約式により同時伐採を回源管理の考え方が出てきた．19 世紀に入るとドイツで避できる．一方 Model II に対しては，下記で対応で保続収穫の考え方が提唱され始め，それに伴い 法正

きる．

林 という保続収穫を可能にする 理想的な 計画モ⎛⎞max(0,t−k)tt⎜ 2 · zij
xist + xjst ≤ 2 · zij≤+ 1, ⎟
⎟
⎜
⎟
⎜
s=0
⎟
⎜
⎟
⎜
,
∀i,
t
=
1,
.
.
.
,
T
∀j
∈
NB
i
⎟
⎜
⎟
⎜ N
T
N
⎟
⎜
B(n)
B(n)
t
⎟
⎜
yij
≤
zij ≤ Kmax ·
yij ,
⎟
⎜
t=1
n=1
⎠
⎝ n=1
∀j ∈ NBi , ∀i

(19)

デルが考案された [31]． 法正林 は，すべての林齢ク
ラスが同じ面積で存在する森林で，理論的には，成長量と等しい材積の収穫が毎期可能になるため，永続的に一定量の木材が供給できる森林であるとしている．
毎回同じ施業を繰り返し，同じ量の収穫を得る森林は必ずしも上記の条件を満たす必要はなく，鈴木氏が提唱した林齢遷移を考慮すれば，その行列の固有値問題
として状態が変わらない， 不動点 を探求することがhModel I の zijに対し，Model II では隣接する i とできる．これを鈴木氏 [1] は 広義の法正林 とした．

j 番目の林分が t 期に伐採されるとき 1 になるバイナ

日本においては，ドイツ林業の影響を強く受け，かつtである．Kmax は計画期間内で実施可リー変数が zijて森林保続思想をもつ研究者，それに対抗する林野行

能な伐採の最大数である．これら (18) あるいは (19)

政の間で 保続か経営か をめぐり長い論争があった

式により集約された林分の同時期でのクラスター伐採

そうである．

などの対応が可能になり，伐採計画の制約，集約の制

20 世紀に入ると，世界大戦の影響もあり，木材需要

約，およびそれらを連結する制約により，最大伐採許が高まり，消費者の動向や木材市場の分析が展開され

容面積制約下での伐採計画最適化が可能になる．なお，た．この時期には，より短期的な経済効率性が管理の

Yoshimoto and Asante [26] は Model I を軸にした指標となり，経済理論に基づいた手法による森林管理

616（50）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ

の計画が検討された [32]．そして，本稿で紹介したよになる．OR 手法の応用にはならないかもしれないが，

うな OR の手法が，実社会のさまざまな森林資源管理

今後の研究の発展に期待したい．

の問題解決に応用され始めた．近年，木材生産から森

ここでは取り上げなかったがほかにも Samuel-

林の環境的な機能に対する社会の関心が変化し，木材

son [34] による空間均衡モデルを用いた世界木材貿易モ

生産とのトレードオフ評価が可能な線形計画法の応用

デル GFPM (Global Forest Products Model) [35] やが進んできた．

日本を中心にした木材貿易モデル JAFSEM (Japanese

世紀の後半 1980 年代に入ると，大規模な皆伐が森Forest Sector Model) [36] も OR 手法を応用したモ林に生息する野生動物の生息地破壊に繋がっているとデルである．今後複雑化する社会ニーズ，自然環境下の批判から，森林の環境的な役割を考慮するうえで，で如何に森林資源を持続的に効率良く管理していくか

伐採地の空間配置を制御することが必要であるというは，OR 研究の成果に依存するところが大きい．北米，考え方が広まった．これに対応すべく混合 0-1 整数計

南米など森林産業が国内の GDP に貢献するところで画法が応用され，隣接する林地の同時期伐採を回避す

は，効率的な資源管理を遂行するためにさまざまなシる制約条件を導入したさまざまな最適化モデルが考案

ステムが開発されコンサルティングも行われている．

された．これらのモデルは，コンピュータ技術の発展

これに対し，ドイツ思想が根強く残る日本においては，と OR 手法の発展に伴い 1990 年代から 2000 年代に

データサイエンスを用いた科学的根拠に依存しない思

かけて，さらに改良・拡張が施され，連続性などの空想的な資源管理がまん延しているように感じる．その間構造も制御できるモデルへと発展していった．そし大きな原因が，日本においては未だエンドユーザーでて，近年，世界各地で頻発する大規模森林火災や病害ある森林所有者，森林組合，行政が容易に利用できる虫の拡散，侵略的外来種の移入・拡散といった空間的ような資源管理システムが存在しないからではと常にな拡散プロセスを伴うさまざまな被害に対して，有効感じている．最適化モデルの開発とともに，そうしたな管理を実施することが重要な課題となってきている．

システムの開発が今後日本における森林資源管理に役

空間拡散プロセス予測モデルと混合 0-1 整数計画法と立つものとなろう．

いった最適化モデルを結合した新たな最適化システムが開発され [33]，管理の空間配置と被害拡散の相互作用を踏まえながら，最適な資源管理を時空間的に探索するような研究が今後ますます展開されていくだろう．
このように，森林資源管理に関する最適化モデルの開発は，森林資源を取り巻く社会ニーズの変化とともに，異なる手法の応用であったり，モデルの改良・拡張を繰り返しながら，発展してきた．
先に取り上げた減反率に関する研究は現在それほど進展がない．その大きな理由が減反率を推定する確率過程によるモデリングにこだわっているからのように思える．実際の伐採の傾向は木材価格，管理費用，利子率といった経済的な要因によって変動するため，単純ポアソン過程から導出されるガンマ分布による減反率の推定は適切ではないことは明らかである．これまで確率過程に基づいた推定法の確立と減反率推定を結び付けた研究が進んできたが，実際の伐採動向を見る
限り，それらは切り離して対応すべき研究内容である．
減反率はその分布を描写することにより伐採率がわかる．仮に実際の伐採率の時系列データと他のファンダ
メンタルズに基づき，パラメトリック，ノンパラメト
リックなどの統計的手法により 減反率分布 を推定

参考文献[1] 鈴木太七，『森林経理学』，朝倉書店，1979.
『木材の生産予測について (I) 』，科学技術庁資[2] 鈴木太七，
源局，1961.
[3] 鈴木太七，『木材の生産予測について (II) 』，科学技術庁
資源局，1963.
[4] 鈴木太七， 遷移確率行列による収穫予定， 日林関東支論，10, pp. 36–38, 1959.
[5] 広嶋卓也， コラム 3.2．減反率と広義の法正林，地域資源の見通し，『森林計画学入門』，田中和博，吉田茂二郎，
白石則彦，松村直人（編），朝倉書店，pp. 43–45, 2020.
[6] 藤掛一郎， 森林資源構成表を用いた減反率の推定， FORMATH, 2, pp. 159–178, 2003.
[7] P. Blandon, “Gentan probability and censored sample theory (I),” Journal of Japanese Forestry Society,73, pp. 187–196, 1991.
[8] A. Yoshimoto, “A new stochastic model for harvesting behavior with application to nonstationary forest
growth and supply,” Canadian Journal of Forest Research, 26, pp. 1967–1972, 1996.
[9] A. Yoshimoto, “Economic analysis of harvesting behavior using the modiﬁed Gentan probability theory,”
Journal of Forest Research, 1, pp. 67–72, 1996.
[10] A. Yoshimoto, “Gentan probability analysis with astate-dependent discrete forest growth model,” Journal of Forest Research, 6, pp. 101–110, 2001.
[11] A. Yoshimoto, “Application of the Logistic, Gompertz, and Richards growth functions to Gentan
probability analysis,” Journal of Forest Research, 6,
pp. 265–272, 2001.

し，分布予測が可能になれば，伐採動向の予測も可能

2022 年 11 月号

Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited. （51）617

[12] J. Martin and P. E. Sendak, Operations Research inForestry: A Bibliography, U.S. Department of Agriculture, Forest Service, Northeastern Forest Experiment
Station, 1973.
[13] K. N. Johnson and H. L. Scheurman, “Techniques
for prescribing optimal timber harvest and investment
under diﬀerent objectives: Discussion and synthesis,”
Forest Science Monograph, 18, 1977.
[14] O. Garcia, “Linear programming and related approaches in forest planning,” New Zealand Journal ofForestry Science, 20, pp. 307–331, 1990.
[15] A. Yoshimoto and P. Asante, “Inter-temporal
aggregation for spatially explicit optimal harvest
scheduling under area restrictions,” Forest Science,
67, pp. 587–606, 2021.
[16] K. N. Johnson, D. B. Jones and B. M. Kent,
A user’s guide to the forest planning model (FORPLAN), Fort Collins, U.S. Department of Agriculture,
Forest Service, Land Management Planning, 1980.
[17] P. A. Arp and D. R. Lavigne, “Planning withgoal programming: A case study for multiple-use offorested land,” Forestry Chronicle, 58, pp. 225–232,
1982.
[18] A. Muray, “Spatial restrictions in harvest scheduling,” Forest Science, 45, pp. 45–52, 2018.
[19] A. Murray and R. L. Church, “Measuring the eﬃcacy of adjacency constraint structure in forest planning models,” Canadian Journal of Forest Research,
25, pp. 1416–1424, 1994.
[20] A. Yoshimoto and M. Konoshima, “Spatially constrained harvest scheduling for multiple harvests by exact formulation with common matrix algebra,” Journal of Forest Research, 21, pp. 15–22, 2016.
[21] M. E. McDill, S. Rebain and J. Braze, “Harvest
scheduling with area-based adjacency constraints,”
Forest Science, 48, pp. 631–642, 2002.
[22] M. Goycoolea, A. T. Murray, F. Barahona, R.
Epstein and A. Weintraub, “Harvest scheduling subject to maximum area restrictions: Exploring exact
approaches,” Operations Research, 53, pp. 490–500,
2005.
[23] M. Constantino, I. Martins and J. G. Borges, “A
new mixed-integer programming model for harvest
scheduling subject to maximum area restrictions,” Operations Research, 56, pp. 542–551, 2008.
[24] 吉本敦，木島真志，栁原宏和， 隣接空間制約による林分
団地化最適パターンの探求， 統計数理，58, pp. 113–126,
2010.
[25] A. Yoshimoto, P. Asante and M. Konoshima, “Op-timal forest stand aggregation and harvest scheduling using compactly formulated integer programming,”
FORMATH, 16, pp. 22–31, 2017.
[26] A. Yoshimoto and A. P. Asante, “A new optimization model for spatially constrained harvest schedulingunder area restrictions through maximum ﬂow problem,” Forest Science, 64, pp. 392–406, 2018.
[27] E. Z. Baskent and G. A. Jordan, “Designing forest management to control spatial structure of landscapes,” Landscape and Urban Planning, 34, pp. 55–
74, 1996.
[28] E. Z. Baskent and S. Keles, “Spatial forest planning: A review,” Ecological Modelling, 188(2-4),
pp. 145–173, 2005.
[29] A. B. Franklin, D. R. Anderson, R. J. Gutiérrez andK. P. Burnham, “Climate, habitat quality, and ﬁtness
in northern spotted owl populations in northwestern
California,” Ecological Monographs, 70, pp. 539–590,
2000.
[30] A. Yoshimoto, and A. P. Asante, “Focal-point
aggregation under area restrictions through spatially constrained optimal harvest scheduling,” Forest
Science, 65, pp. 164–177, 2019.
[31] D. Rose and T. E. Burk, “Development of a model
for simulation of forest regulation techniques,” Agricultural Experiment Station University of Minnesota,
Technical Bulletin, No. 324, 1980.
[32] W. A. Duerr, J. Fedkiw and S. Guttenberg, “Financial maturity: A guide to proﬁtable timber growing,” U.S. Department of Agriculture Technical Bulletin, No. 1146, 1956.
[33] A. Yoshimoto, M. Konoshima, P. Surovy and P.
Asante, “Integer programming approach to control
invasive species spread based on cellular automaton
model,” Natural Resource Modeling, 30, pp. 1–42,
2017.
[34] P. A. Samuelson, “Spatial price equiblium and linear programming,” American Economic Review, 41,
pp. 284–303, 1952.
[35] P. Nepal, J. Buongiorno, C. M. T. Johnson, F.
Prestemon and J. Guo, “Global forest products trademodel,” International Trade in Forest Products: Lumber Trade Disputes, Models and Examples, G. C. vanKooten & L.Voss (eds.), Oxfordshire, pp. 110–141,2021.
[36] A. Yoshimoto and K. Yukutake, “Japanese forest
sector modeling through a partial equilibrium market model,” Journal of Forest Research, 7, pp. 41–48,
2002.

618（52）Copyright © by ORSJ. Unauthorized reproduction of this article is prohibited.

オペレーションズ・リサーチ