The model of stand basal area gross growth on the data of the Estonian Network of Forest Research Plots

Allar Padari; Andres Kiviste; Diana Laarmann; Ahto Kangur

doi:10.2478/fsmu-2023-0007

Sissejuhatus

Kasvava metsa tagavara ja juurdekasvu hindamine on metsateaduses ja -poliitikas aktuaalne teema olnud üle saja aasta. Tänapäeva Eestis on see teema järjest enam hakanud kõnetama erinevate valdkondade esindajaid ja tavainimesi. Terminid, nagu juurdekasv, täisjuurdekasv, tagavara muut, säiliva puistu juurdekasv, surem jne, on kasutuses olnud pikka aega, kuid praeguseks oleme jõudnud olukorda, kus mitme valdkonna esindajad tõlgendavad sageli nende sisu täiesti erinevalt. Eeltoodud loetelus on „juurdekasv“ ja „täis-juurdekasv“ metsanduslikus kontekstis sisult sama tähendusega. Termin „täis-juurdekasv“ on vananenud ja tänapäeval on kasutusel „juurdekasv“ (Vaus, 2005). Täisjuurdekasv või juurdekasv saab olla keskmine, perioodi keskmine või jooksev, millest viimane kirjeldab ühe aasta juurdekasvu. Puistu juurdekasv kirjeldab puistu kõikide puude kasvu, mis olid hinnatava perioodi algul elusad. Kui mõni puu sureb hinnatava perioodi jooksul, siis nende puude mahtusid puistu mahust maha ei lahutata. Kui puistu juurdekasvust lahutada maha surnud puude maht, siis saadakse puistu tagavara muut. Puistu tagavara muut on sama, mis on säiliva puistu juurdekasv ja seda saab arvestada hinnatava puistu kahe ajahetke tagavara vahena – tagavara hinnatava perioodi lõpus miinus tagavara hinnatava perioodi algul. Seega võrdub puistu tagavara või ristlõikepindala juurdekasv puistu tagavara või ristlõikepindala muudu ja mõõtmiste vahepeal surnud (sh raiutud) puude tagavara või ristlõikepindala (puistu suremi) summaga (Gadow et al., 2021; Pretzsch, 2009; Vaus, 2005).

Puistu tagavara juurdekasvu arvutamiseks kasvava metsa takseerkirjelduste and-metel on „Metsa korraldamise juhendi“ (2023) lisas 12 esitatud mudel, mille aluseks on Priit Kohava kolm aastakümmet tagasi koostatud valemid (Kohava, 1992). Kohava koostatud puistu tagavara jooksva juurdekasvu mudel koosneb kahest komponendist: normaalpuistu tagavara juurdekasvu mudelist, mille argumentideks on puistu vanus, boniteet ja puuliik, ning parandus-kordajast K, mis sõltub puistu vanusest, esimese rinde täiusest ja puuliigist. Priit Kohava lähendas normaalpuistu üldtootlikkuse V. Zagrejevi tabelite andmeil (Antanaitis & Zagreev, 1981) Mitscherlichi kasvufunktsiooniga (ForMIS, 2023: kasvuvõrrand 5.03) ning võttis normaalpuistu tagavara juurdekasvuks üldtootlikkuse funktsiooni tuletise. Puistu täiusest tingitud parandus-kordaja tuletas Priit Kohava (Kohava, 1992) Jakase ja Zagrejevi uurimustest. Piisavalt usaldusväärsete Eesti empiiriliste andmete puudumise tõttu pidas Priit Kohava (Kohava, 1992) seda mudelit ajutiseks lahenduseks, tunnistades metsa püsiproovitükkide kordusmõõtmiste vajadust ja puistu juurdekasvu mudeli uuendamist uute mõõtmisandmete saamise järel.

Puistu juurdekasvu mudel peab olema praktikas kasutatav ehk sisendina kasutama metsade traditsioonilisi takseerandmeid (Peuhkurinen et al., 2008). Metsakorralduse mudelite süsteem on pidevas uuendamises (Bravo et al., 2019; Kiviste et al., 2015). Seetõttu on puistu rinnaspindala juurdekasvu valemi eelistamine mahu juurdekasvu ees põhjendatud (Kärki et al., 2000). Rinnaspindala arvestus põhineb diameetrite mõõtmisel, tagavara põhineb lisaks diameetril veel kõrgusel ja mahu või vormiarvu mudelitel.

Puistu rinnaspindala ja tagavara juurdekasvu ei käsitleta metsade modelleerimise uurimistöödes tavaliselt omaette nähtusena, vaid kui ühte komponenti puistu kasvukäiku prognoosivates võrrandisüsteemides. Seejuures teatud puistu kasvu-mudelites juurdekasv ilmutatud valemi kujul polegi esitatud, vaid see tuleb teatud algoritmiga arvutada. Puistu kasvukäigu mudeleid võib sõltuvalt eesmärgist ja kogutud andmete struktuurist koostada erinevatest printsiipidest lähtudes. Metsa majandamise kavandamiseks on enim kasutust leidnud puu kasvu valemitele tuginevad puistu kasvukäigu simulatsioonmudelid (näiteks Soome „MOTTI” (Hynynen et al., 2002) ja Rootsi „Heureka” (Elfving, 2010)), aga endiselt koostatakse ka lihtsamaid, puistu takseertunnuste vahelistele seostele tuginevaid mudeleid (näiteks Allen II et al., 2020; Álvarez-González et al., 2010; García & Ruiz, 2003; Liepa, 2008). Üleminek traditsioonilistelt kasvukäigutabelitelt puu kasvu mudelitele on tekitanud nii optimismi kui ka skeptitsismi (Gadow et al., 2021). On proovitud ka mõlemat tüü-pi mudelite kooskasutust (Cao, 2014; Hevia et al., 2015). Sõltumata mudeli vormist peetakse puistu rinnaspindala juurdekasvu peamisteks mõjufaktoriteks puistu vanust, tihedust ja kasvukoha headust (Sun et al., 2007).

Siinse töö eesmärk oli koostada Eesti enamlevinud puuliikidele üksikpuude mõõtmisandmete alusel uue puistuelemendi juurdekasvu mudel.

Andmed ja metoodika

Puistu juurdekasvu kohta on võimalik saada usaldusväärseid andmeid metsa püsiproovitükkide kordusmõõtmiste abil. Siinses uurimuses on kasutatud mudelite loomiseks Eesti metsa kasvukäigu püsiproovitükkide (KKPRT) võrgustiku (Kiviste et al., 2015) 1995. kuni 2020. aastate mõõtmiste andmeid.

KKPRT võrgustiku andmed katavad enamiku Eesti metsa kasvukohatüüpe ja puistute vanuseklasse, kusjuures proovitükkide asukohad puistutes on valitud juhuslikult. KKPRT proovitükkide raadius on reeglina 15, 20 või 25 m tagamaks vähemalt 100 põhirinde puud proovitükil. Proovitükkidel on mõõdetud kõigi puude asukoha koordinaadid, mis võimaldab kordusmõõtmistel samu puid tuvastada ja uuesti mõõta. Kui proovitükkide mõõtmiste vahelisel ajal on toimunud hooldusraie, siis kordusmõõtmisel proovitüki raadiust suurendatakse, et proovitükil oleks 100 puu nõue täidetud.

Proovitükkidel mõõdeti iga puu rinnasdiameeter kahes suunas ja tuvastati puu seisund (rinne, elus, surnud, raiutud, vigastused). Mudelpuudel mõõdeti ka puu kõrgus ja võra alguse kõrgus. Puistu vanuseks proovitükil võeti I rinde vanus metsakorralduse andmeist, mida kontrolliti ja vajadusel korrigeeriti ning täiendati puursüdamike aastarõngaste lugemise teel.

Puistu takseertunnused (rinnaspindala, täius, boniteet) arvutati proovitüki mõõtmisandmete põhjal „Metsa korraldamise juhendi” (2023) järgi, kusjuures täpselt proovitüki piirile jäävate puude rinnaspindalast võeti arvesse pool. Mõõtmiste vahelisel perioodil eluspuude, raiutud ja surnud puude juurdekasv arvestati eraldi puuliikide ja rinnete osas.

Täius ja boniteet H₁₀₀ arvutati „Metsa korraldamise juhendi“ (2009) lisades 11 ja 10 esitatud valemite alusel. Ristlõikepindala arvutamiseks kasutati rinnasdiameetrite mõõtmisandmeid.

Kuna proovitükkide kordusmõõtmised ei pruugi toimuda esimese mõõtmisega samal kuupäeval, siis arvestati puistu vanuse arvutamisel ka mõõtmise kuupäeva. Puistu vanusele kasvu sesoonsuse arvutamiseks kasutati mudeleid, mille koostamine on kirjeldatud Kiviste et al. (2022) artiklis. Mudeli kuju määramiseks kasutati Vitase (2011) andmeid. Kuna Eesti asub Leedust põhja pool, siis eeldati, et puude radiaalkasv algab Eestis hiljem ja lõpeb varem kui Leedus. Puude radiaalkasvu alguse ja lõpu määramiseks Eesti laiuskraadile kasutati lisaks Leedu (Vitas, 2011) andmetele Soome (Mäkinen et al., 2008) vastavaid uurimistulemusi (Kiviste et al., 2022).

Kuna eesmärk oli arvutada puistu ristlõikepindala juurdekasvu, siis osutus oluliseks kasutada kõikide puude andmeid, mis olid mõõtmisperioodi algul elusad. Puistu juurdekasvu arvutamisel on oluline, millal puu suri või kui palju puu jõudis enne suremist kasvada. Antud töös eeldati, et perioodi jooksul oli puude suremus ajaliselt ühtlane ning puid kasvatati edasi kuni hinnatava perioodi keskele. Surnud puude diameetri juurdekasvu arvutamiseks kasutati metoodikat, mille järgi peenemate puude juurdekasv oli väiksem ja jämedate puude juurdekasv suurem. Selleks loodi metoodika, kus algul leitakse kõikidele esimese rinde elusatele puudele diameetri kasvu kirjeldav valem, mis läbib täpselt looduses mõõdetud mõlemat diameetrit:

1

d = a \cdot {(\frac{A}{A + 1})}^{b},

kus d on puu diameeter (cm), A on puu vanus (aasta) ning a ja b on valemi parameetrid. Igale puule proovitükil leitakse valemi parameetrid a ja b sirgeks teisendatud võrrandit kasutades:

2

\ln d = \ln a + b \cdot \ln (\frac{A}{A + 1}),

millele omakorda parameetrite b ja a arvutamisel kasutatakse perioodi alguses mõõdetud diameetrit (d₁) ja hinnatud vanust (A₁) ning perioodi lõpus mõõdetud diameetrit (d₂) ja hinnatud vanust (A₂):

3

b = \frac{\ln (d_{2}) - \ln (d_{1})}{\ln (\frac{A_{2}}{A_{2} + 1}) - \ln (\frac{A_{1}}{A_{1} + 1})},

4

a = e^{\ln (d_{1}) - b \cdot \ln (\frac{A_{1}}{A_{1} + 1})} .

Pärast elusatele puudele parameetrite a ja b arvutamist teostati regressioonanalüüs valemile, mille kasutamisel sai arvutada puu diameetri suvalises vanuses. Kasutades regressioonanalüüsi, leiti parameetri a sõltuvus puu diameetrist astmefunktsiooniga ning parameeter b sõltuvus diameetrist leiti logaritm-funktsiooniga:

5

\begin{matrix} a = a_{0} \cdot d^{a_{1}} \\ (lineaarne kuju : \ln a = ln a_{0} + a_{1} \cdot \ln d), \end{matrix}

6

b = b_{0} + b_{1} \cdot \ln d,

kus a ja b on valemi 1 parameetrid, a₀, a₁, b₀ ja b₁ on valemite 5 ja 6 parameetrid ning d on puu diameeter.

Sellisel viisil leiti parameetrid a₀, a₁, b₀ ja b₁ igale proovitükile ja igale esimese rinde puistuelemendile. Surnud puule arvutati seejärel valemitega 5 ja 6 parameetrid a ja b, kusjuures diameetrina kasutati perioodi algul mõõdetud diameetrit. Surnud puudele arvestati kasvamiseks pool perioodi ning kasutades valemis 1 poolel kasvuperioodil olnud puu vanust ning parameetreid a ja b, leiti surnud puude viimane diameeter.

Puistu ristlõikepindala juurdekasvu arvutamiseks eristati ristlõikepindalad rinnete kaupa. Edasisel modelleerimisel kasutati esmajoones esimese rinde ristlõikepindalasid. Puistu juurdekasvu mudeli parandi arvutamisel kasutati ka koefitsienti, mis näitab, kui suure osa moodustab esimese rinde ristlõikepindala kõikide eluspuude rinnete ristlõikepindalade summast. Puistu juurdekasvu põhivalemi arvutamise aluseks võeti esimese rinde ristlõikepindala. Kuna puistu ristlõikepindala juurdekasv on reeglina kasvav ja seda kirjeldav funktsioon ei ole lineaarne, siis järelikult on puistu ristlõikepindala juurdekasv keskmiselt muutuv ning esimese aasta juurdekasv ei ole sama, mis on perioodi keskmine ristlõikepindala juurdekasv. Seega kasutati esimese aasta juurdekasvu arvutamiseks valemit, mille parameetrite leidmiseks kasutati regressioonanalüüsi:

7

G_{1} = c_{0} \cdot {(\frac{A}{A + 1})}^{c_{1}},

kus G₁ on esimese rinde ristlõikepindala vanusel A ning c₀ ja c₁ on valemi parameetrid. Valemi kuju valiti selline, mis algaks koordinaatide alguspunktist (0,0). Valemi konstantide leidmiseks kasutati regressioonanalüüsi, mille sisendiks olid perioodi alguse ja lõpu mõõtmise andmeid. Parameetrite c₀ ja c₁ hinnangud leiti analoogselt nagu valemile 1 leiti hinnangud parameetritele a ja b. Perioodi alguse ja lõpu vahe oli enamasti neli kuni kuus aastat. Puistu esimese aasta ristlõikepindala juurdekasv arvutati järgmiselt:

8

z G = a_{0} \cdot {(\frac{A_{1} + 1}{A_{1} + 2})}^{a_{1}} - G_{1},

Puistu juurdekasvu hindamisel kasutati ka esimese rinde ristlõikepindala suhet kõikide elusate puudega rinnete ristlõikepindalasse ehk G₁/G_sum. Puistu juurdekasvu mudelite koostamisel valiti välja puistud, mille kõrgus oli vähemalt 6 m ning mille ristlõikepindalade suhe G₁/G_sum oli vähemalt 0,95. Samuti jäid välja need mõõtmispaarid, millel jäi kahe mõõtmise vahele raie, kui väljaraie oli suurem kui 5%. Hariliku männi (Pinus Sylvestris L., edaspidi: mänd), hariliku kuuse (Picea abies (L.) H. Karst., edaspidi: kuusk) ja halli lepa (Alnus incana (L.) Moench) mudelite koostamisel kasutati vaid proovitükke, kus vastav puuliik oli ka enamuspuuliik. Kase (eristamata: arukask (Betula pendula Roth) ja sookask (Betula pubescens Ehrh.), edaspidi: kask), hariliku haava (Populus tremula L., edaspidi: haab) ja sanglepa (Alnus glutinosa (L.) Gaertn.) puhul võeti valemi koostamisel kasutusse ka kõrvalpuuliigid. Puuliigi osakaalu miinimumiks oli 10%. Kõrvalpuuliigi andmeid kasutati nende puuliikide puhul seetõttu, et enamuspuuliigina ei oleks analüüsiks olnud piisavalt andmeid. Ka halli lepa puhul oli liiga vähe andmeid, kuid halli lepa puistuelementide juurdekasvud enamuspuuliigina ja kõrvalpuuliigina osutusid liiga erinevaks.

Kasutusele võetud andmete esinemissagedused boniteediklasside ja täiusklasside kaupa on toodud tabelites 1 kuni 6 vastavalt puuliikidele mänd, kuusk, kask, haab, sanglepp ja hall lepp.

Enne modelleerimist jagati esimese rinde puistuelementide ristlõikepindalade juurdekasvud puistuelemendi ristlõikepindala osakaaluga kogu esimese rinde ristlõikepindalast. Puistu ristlõikepindala juurdekasvu ja vanuse vahelise seose leidmiseks tehti regressioonanalüüs. Mudeli kuju valiti selline, mille alguspunkt oleks (0,0) ning vanuse suurenedes oleks algul funktsioon tõusev, hiljem langev ning funktsioon ei annaks vanuse positiivse väärtuse korral ühelgi juhul negatiivset väärtust. Nendele omadustele vastav ja selles töös sobilik valemi kuju on järgmine:

9

z G = e^{b_{0} + b_{1} \cdot \ln A + b_{2} \cdot {(\ln A)}^{2}},

mis lineaarsele kujule teisendatult on järgmine:

10

\ln z G = b_{0} + b_{1} \cdot \ln A + b_{2} \cdot {(\ln A)}^{2} .

Tabel 1.

Männi enamusega proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of Scots pine dominated plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative density, %	Boniteediklass
	Site index class
	Ib	Ia	I	II	III	IV	V	Va	Vb
10			3 (51–55)
30		2 (76–81)		2 (61–66)			2 (35–38)
50	1 (62)	4 (34–53)	26 (18–126)	23 (20–106)	21 (25–120)	8 (29–101)	7 (75–191)	1 (62)
70	2 (21–44)	22 (21–86)	102 (26–103)	87 (23–131)	30 (23–146)	39 (35–126)	30 (61–161)	4 (67–85)	1 (76)
90		39 (21–81)	75 (18–130)	78 (20–106)	41 (26-116)	37 (26-241)	27 (66–246)
110		8 (35–89)	24 (24–96)	28 (24–110)	12 (21–111)	4 (26–111)	1 (176)
130			5 (20–38)	5 (45–95)	3 (26–50)	3 (111–121)
150			1 (24)
Kokku Total	3 (21–62)	76 (21–89)	236 (18–130)	223 (20–131)	107 (21–146)	91 (26–241)	67 (35–246)	5 (62–85)	1 (76)

Tabel 2.

Kuuse enamusega proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of Norway spruce dominated plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative density, %	Boniteediklass
	Site index class
	Ib	Ia	I	II	III	IV
30		3 (19–28)	1 (21)		2 (26)
50	1 (34)	4 (24–60)	3 (21–36)	8 (25–92)	1 (66)	1 (63)
70	10 (15–50)	24 (26–66)	21 (26–66)	4 (31–82)	3 (31–91)	2 (68–73)
90	11 (19–50)	51 (25–80)	32 (31–89)	8 (35–61)	1 (66)
110	13 (20–43)	43 (26–56)	23 (30–80)	12 (35–105)	1 (51)	1 (101)
130	4 (17–40)	6 (34–50)	11 (31–51)	3 (91–91)	3 (41–51)
150		2 (31–33)	1 (40)
Kokku Total	39 (15–50)	133 (19–80)	92 (21–89)	35 (25–105)	11 (26–91)	4 (63–101)

Tabel 3.

Kase (osakaal ≥ 10%) proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of birch (share ≥ 10%) plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative density, %	Boniteediklass
	Site index class
	Ib	Ia	I	II	III	IV	V – Vb
10		1 (51)	2 (13)
30		1 (13)	2 (13)		1 (26)	2 (23–28)	2 (21–38)
50		1 (13)	10 (13–54)	3 (25–31)	6 (23–36)	2 (23–31)
70	4 (21–40)	22 (13–76)	36 (13–95)	15 (23–71)	14 (21–63)	4 (36–73)	3 (46–51)
90		37 (13–101)	29 (15–76)	41 (21–74)	15 (21–73)	5 (28–53)	1 (41)
110	3 (28–41)	12 (17–51)	19 (20–86)	31 (21–105)	11 (25–51)	4 (36–58)
130	2 (17–22)		6 (16–46)	10 (28–95)	6 (30–51)
150		1 (31)	1 (24)	2 (33–41)	4 (30–51)
170				1 (45)
190					1 (35)
Kokku Total	9 (17–41)	75 (13–101)	105 (13–95)	103 (21–105)	58 (21–73)	17 (23–73)	6 (21–51)

Tabel 4.

Haava (osakaal ≥ 10%) proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of aspen (share ≥ 10%) plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative density, %	Boniteediklass
	Site index class
	Ib	Ia	I	II	III	IV	V – Vb
30			3 (31–41)
50	4 (55–81)	12 (33–136)	11 (31–80)			1 (31)	2 (96–101)
70	10 (31–101)	44 (25–101)	50 (32–91)	11 (23–106)	3 (27–96)	2 (36–90)	1 (49)
90	2 (70–101)	40 (26–92)	24 (36–80)	6 (41–111)	2 (32–40)	3 (35–53)
110	2 (60–66)	6 (35–86)	12 (30–100)	8 (27–105)	2 (41–114)	1 (58)
130		1 (47)	5 (16–56)	2 (56–61)			2 (80–111)
150				1 (65)
Kokku Total	18 (31–101)	103 (25–136)	105 (16–100)	28 (23–111)	7 (27–114)	7 (31–90)	5 (49–111)

Tabel 5.

Sanglepa (osakaal ≥ 10%) proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of common alder (share ≥ 10%) plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative density, %	Boniteediklass
	Site index class
	Ib	Ia	I	II	III	IV	V – Vb
10
30							1 (21)
50			4 (25–50)	1 (25)	4 (23–56)	1 (63)
70		3 (40–52)	6 (30–66)	5 (28–70)	1 (21)	2 (68–73)
90		2 (21–37)	12 (15–60)	15 (21–56)	2 (21–40)	4 (28–78)
110	1 (37)	8 (26–46)	16 (20–60)	10 (27–65)	8 (25–56)	1 (58)
130			7 (30–56)	6 (28–71)	2 (30)
150				4 (33–65)
170				1 (45)	1 (30)
190					1 (35)
Kokku Total	1 (37)	13 (21–52)	45 (15–66)	42 (21–71)	19 (21–56)	8 (28–78)	1 (21)

Tabel 6.

Halli lepa proovitükkide mõõtmispaaride esinemissagedused boniteediklasside ja täiusklasside kaupa (sulgudes on esimese mõõtmise vanuse vahemik).

Frequencies of measurement pairs of grey alder (share ≥ 10%) plots by site index classes and relative density classes (in parentheses is the age range of the first measurement).

Täius Relative %	Boniteediklass
	Site index class
	I	II	III
70		2 (21–32)	1 (27)
90	1 (15)	1 (25)	3 (23–37)
110	1 (46)	4 (26–42)
130		1 (40)	1 (40)
150			3 (26–30)
170		1 (30)
Kokku Total	2 (15–46)	9 (21–42)	8 (23–40)

Selle valemikuju omadused:

Algab punktist 0,0;
Maksimumpunkti asukohta x-telje suhtes (vanus) kirjeldab parameetrite b₁ ja b₂ hinnangute jagatis;
Maksimaalne juurdekasvu väärtus sõltub parameetrite b₀ ja b₁. hinnangute jagatisest;
Juurdekasvu kahanemiskiirus pärast ekstreemumpunkti sõltub parameetri b₁ hinnangu väärtusest;
Vanuse lähenemisel lõpmatusele läheneb juurdekasv 0-le.

Regressioonanalüüs tehti puuliikide ja boniteediklasside kaupa. Seejärel analüüsiti valemi parameetrite hinnangute sõltuvust boniteedist. Kuna valemi 9 ekstreemumpunkt sõltub avaldisest b₁/b₂, siis esmalt analüüsiti b₁/b₂ sõltuvust boniteedist H₁₀₀:

11

\frac{b_{1}}{b_{2}} = k_{0} + k_{1} \cdot H_{100},

Edasiseks analüüsiks arvutati puuliikide ja vanuseklasside kaupa puistu ristlõikepindala juurdekasvu sõltuvus boniteedist H₁₀₀ valemiga:

12

z G = m \cdot H_{100},

kusjuures parameeter m leiti valemiga:

13

m = \frac{\sum \frac{z G}{H_{100}}}{n}

Kuna igale vanuseklassile leiti eraldi parameetri m hinnang, siis tulemuseks saadi punktiparv, mis võeti kasutusse modelleerimise järgmistes etappides. Esmalt võeti kasutusele valem 11, mille kasutuselevõtuga sai valemi 10 konstantide arvu vähendada kahele:

14

\ln z G = b_{0} + b_{1} \cdot (\ln A + \frac{{(\ln A)}^{2}}{k_{0} + k_{1} \cdot H_{100}}),

Leides lineaarse regressiooniga valemi 14 parameetrite b₀ ja b₁ hinnangud kõikidele boniteediklassidele, mille järel leiti valemi parameetri b₀ seos H₁₀₀-st (erinevatel puuliikidel osutus parimaks erinev valemikuju):

15

b_{0, M A, K U} = o_{0} + o_{1} \cdot H_{100} + o_{0} - H_{100}^{2},

16

b_{0, K S, H B, L M, L V} = o_{0} + o_{1} - \ln (H_{100}) .

Teades parameetrite b₁ ja b₂ jagatise seost boniteediga H₁₀₀ ja parameetri b₀ seost boniteediga H₁₀₀, on vaja leida veel parameetri b₁ seos boniteediga H₁₀₀. Selleks kasutati valemite teisendusi ning saadi üheparameetriline valem:

17

\ln z G - b_{0} = b_{1} \cdot (\ln A + \frac{{(\ln A)}^{2}}{k_{0} + k_{1} \cdot H_{100}}),

kus parameeter b₀ arvutamiseks kasutatakse olenevalt puuliigist valemit 15 või 16. Leidnud igale boniteediklassile valemi 17 parameetrile b₁ hinnangud saadi regressioonanalüüsi kasutades puuliikide kaupa parameetri b₁ seos boniteediga H₁₀₀:

18

b_{1, M A, K U} = r_{0} + r_{1} \cdot H_{100} + r_{2} \cdot H_{100}^{2};

19

b_{1, K S, H B} = r_{0} + r_{1} \cdot \ln (H_{100});

20

b_{1, L M, L V} = r_{0} .

Pärast põhivalemi 9 parameetrite hindamist hinnati iga puuliigi valemi parandit K (tegeliku ristlõikepindala ja valemiga 9 arvutatud puistu ristlõikepindala juurdekasvu jagatis). Parandi vajadus tuleneb sellest, et puistute tihedus (täius) on erinev ning eespool kirjeldatud mudeli koostamise metoodikas ei ole arvestatud puistu täiusega. Kuna algandmetel ei pruugi kõikides vanuseklassides ja kõikides boniteediklassides olla täius ühtlase proportsiooniga, siis parandi K modelleerimisel kontrolliti puistu vanuse ja boniteedi kaasamise vajadust. Sellest tulenevalt osutus valemikuju põhivalemi parandi K modelleerimisel puuliigiti mõnevõrra erinevaks:

21

K_{M A} = {(H_{100})}^{S_{1}} \cdot {(\frac{A}{A + 1})}^{S_{2}} \cdot {(\frac{T}{T + 10})}^{S_{3}} \cdot {(\frac{G_{1}}{G_{sum}})}^{S_{4}},

22

K_{K U, L M} = {(\frac{H_{100}}{H_{100} + 10})}^{S_{1}} \cdot {(\frac{A}{A + 1})}^{S_{2}} \cdot {(\frac{T / 100}{T / 100 + 10})}^{S_{3}} \cdot {(\frac{G_{1}}{G_{s u m}})}^{S_{4}},

23

K_{K S, H B} = {(\frac{H_{100}}{H_{100} + 10})}^{S_{1}} \cdot {(\frac{A}{A + 1})}^{S_{2}} \cdot {(\frac{T}{T + 10})}^{S_{3}} \cdot {(\frac{G_{1}}{G_{s u m}})}^{s_{4}},

24

K_{L V} = s_{3} \cdot T .

kus K – põhivalemi parand; H₁₀₀ – puistu boniteet, m; A – puistu vanus, aasta; T – puistu esimese rinde täius, %; G₁ – esimese rinde rinnaspindala, m²/ha; G_sum – kõikide eluspuudega rinnete summaarne rinnaspindala, m²/ha ning s₁, s₂, s₃ ja s₄ – valemi parameetrid.

Männi puhul teostati regressioonanalüüs kahes etapis. Esimeses etapis oli valimis vaid proovitükid, mille I rinde ristlõikepindala kõikide elusate puude ristlõikepindalast oli vähemalt 0,95. Teises etapis olid valimis kõik proovitükid. Kahe etapi regressioonimudelid lineaarsele kujule teisendatuna on järgmised:

25

\ln (\frac{z G_{t e g e l i k}}{z G_{v a l e m i g a}}) = s_{1} \cdot \ln (H_{100}) + s_{2} \cdot \ln (\frac{A}{A + 1}),

26

\ln (\frac{z G_{t e g e l i k}}{z G_{v a l e m i g a} \cdot H_{100}^{s_{1}} \cdot {(\frac{A}{A + 1})}^{s_{2}}}) = s_{3} \cdot \ln (\frac{T}{T + 10}) + s_{4} \cdot \ln (\frac{G_{1}}{G_{s u m}}) .

Tulemused ja arutelu

Puistu ristlõikepindala juurdekasvu mudeli parameetrid ja headus

Puuliikide ja boniteediklasside kaupa hinnatud valemi 9 parameetri hinnangud on toodud mändidele, kuuskedele, kaskedele ja haabadele tabelites 7–10.

Tabel 7.

Männi enamuspuuliigiga proovitükkide ristlõikepindala juurdekasvu valemi parameetrid boniteediklasside ja täiusklasside kaupa.

Boniteediklass	Täius, %	Keskmine H₁₀₀, m	b₀	b₁	b₂	R²	p
Ib	60–80	36,35	-4,8644	4,2620	-0,8167	1,000	–
Ia	60–80	33,05	-4,0588	2,9231	-0,5012	0,865	<0,0001
	80–100		-3,8472	2,5423	-0,4192	0,742	<0,0001
	100–120		-3,4670	1,8774	-0,2762	0,961	<0,0001
I	40–60	29,39	-3,7465	2,4900	-0,4225	0,857	<0,0001
	60–80		-4,1953	3,0514	-0,5281	0,564	<0,0001
	80–100		-3,8412	2,6496	-0,4446	0,723	<0,0001
	100–120		-4,2458	3,1550	-0,5448	0,350	0,0087
	120–140		-4,6025	3,7939	-0,6914	0,908	0,0281
II	40–60	25,76	-3,8764	2,5730	-0,4277	0,853	<0,0001
	60–80		-3,9056	2,6448	-0,4419	0,778	<0,0001
	80–100		-3,8961	2,6447	-0,4389	0,691	<0,0001
	100–120		-3,9584	2,7186	-0,4446	0,854	<0,0001
III	40–60	21,84	-3,6102	2,1957	-0,3505	0,771	<0,0001
	60–80		-3,6258	2,3200	-0,3826	0,777	<0,0001
	80–100		-3,8751	2,6338	-0,4412	0,757	<0,0001
	100–120		-4,1997	3,1175	-0,5289	0,849	<0,0001
	120–140		-4,6315	3,8544	-0,7136	1,000	0,0057
IV	40–60	17,40	-4,00455	2,78573	-0,47857	0,943	0,0002
	60–80		-3,47127	1,95034	-0,30184	0,649	<0,0001
	80–100		-3,30865	1,85014	-0,28675	0,500	<0,0001
	100–120		-4,08709	2,90979	-0,48152	0,997	0,0028
V	40–60	14,36	-3,58913	2,08848	-0,33732	0,937	0,0010
	60–80		-3,41914	1,80899	-0,2765	0,555	<0,0001
	80–100		-3,72392	2,29957	-0,36709	0,824	<0,0001

Tabel 8.

Kuuse enamuspuuliigiga proovitükkide juurdekasvu valemi parameetrid boniteediklasside kaupa.

Boniteediklass	Keskmine H₁₀₀, m	b₀	b₁	b₂	R²	p
Ib	37,10	-4,5447	3,6923	-0,6807	0,869	<0,0001
Ia	33,37	-4,0200	2,8075	-0,4702	0,797	<0,0001
I	29,94	-3,9464	2,6806	-0,4379	0,857	<0,0001
II	25,23	-3,9778	2,7453	-0,4503	0,809	<0,0001
III	21,25	-3,7942	2,4168	-0,3726	0,916	<0,0001

1 Kuna IV boniteedi kuusikuid oli vaid 4, siis neile regressioonanalüüsi ei teostatud (parameetrite p-väärtused olid > 0,05).

Tabel 9.

Kase puistuelementide juurdekasvu valemi parameetrid boniteediklasside kaupa.

Boniteediklass	Keskmine H₁₀₀, m	b₀	b₁	b₂	R²	p
Ib	37,04	-2,0000	1,8882	-0,3659	0,369	0,2517
Ia	33,13	-2,1401	2,1330	-0,4442	0,494	<0,0001
I	29,71	-3,3392	2,6574	-0,4934	0,562	<0,0001
II	25,45	-2,8802	2,0876	-0,3734	0,363	<0,0001
III	21,83	-2,4755	2,1359	-0,4222	0,258	0,0003
IV	17,84	-13,104	7,6270	-1,1250	0,138	0,3543

Tabel 10.

Haava puistuelementide juurdekasvu valemi parameetrid boniteediklasside kaupa.

Boniteediklass	Keskmine H₁₀₀, m	b₀	b₁	b₂	R²	p
Ib	36,98	-1,4583	1,6081	-0,3280	0,928	<0,0001
Ia	33,22	-2,2114	1,4658	-0,2463	0,834	<0,0001
I	29,96	-2,0989	1,4895	-0,2601	0,714	<0,0001
II	25,42	-1,9248	1,5237	-0,2761	0,933	<0,0001
III	22,56	-1,3016	1,6505	-0,3393	0,968	0,0002
IV	16,90	-1,1009	1,6678	-0,3712	0,993	<0,0001
V	14,31	-3,5342	1,1920	-0,1184	0,987	0,0132

Tabelites 7–10 avaldatud parameetrite hinnangute b₁ ja b₂ jagatisele leiti seos boniteedist H₁₀₀ valemiga 11. Regressioonanalüüsi tulemusena saadi valemid, mis on kujutatud joonisel 1 (kõikide valemite puhul oli p < 0,05). Sanglepa ja halli lepa puhul avaldise b₁/b₂ usaldusväärset seost boniteediga ei leitud, sest selleks oli andmeid liiga vähe. Samas konstandina saadi b₁ ja b₂ jagatiseks sanglepal -4,9644 ja hallil lepal -6,4998.

Valemi 5 parameetrite suhte (b₂/b₃) sõltuvus puistu boniteedist ja puuliigist.

Valemiga 12 arvutati igale puuliigile kõikidele esinenud vanuseklasside kaupa puistu ristlõikepindala juurdekasvu lineaarne seos boniteediga. Tulemused on esitatud tabelis 11 ning joonisel 2. Halli lepa puhul oli kasutada vaid 19 mõõtmispaari andmeid ning seetõttu arvutati valemile 11 parameetri a hinnangud kõigi esinenud boniteedi H₁₀₀ kohta. Halli lepa tulemused on näha joonisel 3.

Männikute (A), kuusikute (B), kaasikute (C), haavikute (D) ja sanglepikute (E) ristlõikepindalade juurdekasv erinevates vanustes boniteediklasside kaupa.

Tabel 11.

Puistu rinnaspindala juurdekasvu sõltuvus boniteedist vanuseklasside kaupa.

Vanuseklass	Mänd		Kuusk		Kask		Haab		Sanglepp
Vanuseklass	n	m	n	m	n	m	n	m	n	m
5	1	0,01355
15	17	0,05550	6	0,03902	6	0,03967	1	0,03548	1	0,04361
25	74	0,05016	33	0,03703	66	0,03788	6	0,04758	23	0,03377
35	78	0,03750	96	0,03513	116	0,03648	22	0,04238	33	0,02882
45	101	0,02959	105	0,02980	110	0,02861	39	0,03176	17	0,02572
55	147	0,02336	32	0,02530	104	0,01724	47	0,02790	24	0,02426
65	109	0,02306	16	0,02098	65	0,01416	39	0,02134	13	0,01867
75	117	0,02182	9	0,02049	41	0,01160	38	0,01983	3	0,01829
85	70	0,02265	10	0,02218	36	0,01093	31	0,01827
95	45	0,02051	4	0,02300	20	0,01248	13	0,02136
105	35	0,01754	4	0,02155	20	0,01055	11	0,01977
115	17	0,01648			7	0,00378	3	0,02806
125	9	0,02063
135	3	0,01031					1	0,01158
145	3	0,02290
155	4	0,01667
165	4	0,02248
175	3	0,01218
185	3	0,01154
195	1	0,01153
235	2	0,00604
245	2	0,01161

Hall-lepikute ristlõikepindalade juurdekasv erinevates vanustes boniteediklasside kaupa.

Algandmetena joonisel 2 kujutatud andmeid kasutades leiti valemi 14 parameetritele b₀ ja b₂ regressioonanalüüsi kasutades hinnangud. Regressioonanalüüsil saadud valemid koos graafiliste kujutistega on mändidele, kuuskedele, kaskedele, haabadele, sangleppadele ja hall leppadele vastavalt joonistel 4–9. Joonistelt 4–9 saadud valemite vabaliikmete b₀ hinnangute seosed boniteediga H₁₀₀ on mändidele, kuuskedele, kaskedele, haabadele, sangleppadele ja hall leppadele esitatud valemile vastavalt joonistel 10–15. Saadud b₀ valemeid kasutatakse b₀ arvutamiseks põhivalemi (9) kasutamisel. Selleks et valemile 9 saada kõikide parameetrite hinnangud, siis on vaja veel modelleerida parameetri b₁ sõltuvus boniteedist H₁₀₀. Selleks kasutati regressioonanalüüsi ja valemi 17 kuju. Mändidele, kuuskedele, kaskedele ja haabadele leitud valemi 17 parameetri hinnangud on vasta-valt joonistel 16–19. Sanglepa ja halli lepa puhul parameetri b₁ seost boniteediga H₁₀₀ ei olnud ja seega kasutatakse nende puuliikide puhul b₁ puhul väärtusi vastavalt 1,246 ja 11,9181.

Männikute ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Kuusikute ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Kase puistuelementide ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Haava puistuelementide ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Sanglepa puistuelementide ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Halli lepa puistuelementide ristlõikepindalade juurdekasvu logaritmi sõltuvus puistu vanuse ja boniteedi avaldisest (valemi 14) boniteediklasside kaupa.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist männikute puhul.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist kuusikute puhul.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist kase puistuelemendi puhul.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist haava puistuelemendi puhul.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist sanglepa puistuelemendi puhul.

Valemi 9 parameetri b₀ sõltuvus puistu boniteedist halli lepa puistuelemendi puhul.

Männikute ristlõikepindalade juurdekasvu avaldise sõltuvus vanuse ja boniteediklassi avaldisest (valem 17).

Kuusikute ristlõikepindalade juurdekasvu avaldise sõltuvus vanuse ja boniteediklassi avaldisest (valem 17).

Kase puistuelementide ristlõikepindalade juurdekasvu avaldise sõltuvus vanuse ja boniteediklassi avaldisest (valem 17).

Haava puistuelementide ristlõikepindalade juurdekasvu avaldise sõltuvus vanuse ja boniteediklassi avaldisest (valem 17).

Joonistel 20–23 on toodud joonistel 16–19 kujutatud regressioonivalemite parameetrite hinnangute seosed boniteedist H₁₀₀.

Valemi 9 parameetri b₁ sõltuvus puistu boniteedist männikute puhul.

Valemi 9 parameetri b₁ sõltuvus puistu boniteedist kuusikute puhul.

Valemi 9 parameetri b₁ sõltuvus puistu boniteedist kase puistuelemendi puhul.

Valemi 9 parameetri b1 sõltuvus puistu boniteedist haava puistuelemendi puhul.

Seejärel modelleeriti põhivalemi parand K, mille valemikuju sõltub puuliigist ja on esitatud metoodika peatükis valemite 21–24 all. Regressioonanalüüsi tulemused on toodud tabelites 12 ja 13. Kuna männi puhul teostati parandi K parameetrite hinnangute leidmiseks regressioonanalüüs kahes etapis, siis tabelis 12 on toodud ka determinatsioonikordaja ja p-väärtus kaks korda (valemid 25 ja 26). Mudeli kasutamisel arvutatakse puistu ristlõikepindala juurdekasv, kasutades valemit 9, mille tulemus korrutatakse juurdekasvu parandiga K (sõltuvalt puuliigist valemid 21–24).

Tabel 12.

Puistu juurdekasvuvalemi K parameetrite hinnangud kuuse, männi ja kase puistuelementidele.

Parameeter	Kuusk (valem 22)		Mänd (valem 21)		Kask (valem 23)
Parameeter	hinnang	p	hinnang	p	hinnang	p
s₁	-3,0160	<0,0001	-0,1216	<0,0001	-0,9095	0,0301
s₂	-7,3726	<0,0001	-11,5142	<0,0001	-2,7151	0,5638
N			810
R²			0,2042
p			<0,0001
s₃	0,4279	<0,0001	0,6164	<0,0001	2,6006	<0,0001
s₄	0,4977	<0,0001	0,7270	<0,0001	1,9882	<0,0001
N	495		1173		835
R²	0,4039		0,0852		0,1167
p	<0,0001		<0,0001		<0,0001

Tabel 13.

Puistu juurdekasvuvalemi K parameetrite hinnangud haava, sanglepa ning halli lepa puistuelementidele.

Parameeter	Haab (valem 23)		Sanglepp (valem 22)		Hall lepp (valem 24)
Parameeter	hinnang	p	hinnang	p	hinnang	p
s₁	-2,1190	<0,0001	-2,6860	<0,0001
s₂	0		-4,8688	0,2460
s₃	5,9814	<0,0001	0,4772	<0,0001	0,0116	<0,0001
s₄	0		0,0511	0,9490
N	274		130		19
R²	0,3414		0,2530		0,9287
p	<0,0001		<0,0001		<0,0001

Valemi vigade analüüsi (mõõdetud puistu juurdekasv miinus mudeliga arvutatud puistu juurdekasv) visuaalsed tulemused on toodud järgmistes lisades: lisas 1–3 on männi, 4–6 kuuse, 7–9 kase, 10–12 haava, 13–15 sanglepa ning 16–18 halli lepa valemi vead võrrelduna vastavalt boniteedist H₁₀₀, vanusest ja täiusest.

Olemasoleva puistu juurdekasvu mudeli ja koostatud puistu juurdekasvu mudeli võrdlus

Siinse töö eesmärk oli koostada puistu juurdekasvu mudel kasvukäigu püsiproovitükkide kordusmõõtmistel kogutud andmete alusel. Kuna praegune Eesti tingimustes kasutamiseks olev puistu juurdekasvu mudel („Metsa korraldamise juhend“, 2023: lisa 12) ei ole koostatud Eesti püsiproovitükkide andmeil, siis oli põhjust arvata, et mudeliga arvutatud tulemused võiksid anda tegelikkusega võrreldes liiga suure vea. Joonisel 24 on esitatud puuliikide kaupa puistu tagavara juurdekasvu arvutamistulemuste erinevused (selles töös koostatud mudeliga arvutatud puistu tagavara juurdekasv miinus praegu kehtiv puistu tagavara juurdekasv). Joonisel 24 on violetse joonega toodud puistu juurdekasvude erinevused V boniteediklassi (H₁₀₀ = 13,5 m), sinise joonega III boniteediklassi (H₁₀₀ = 21,5 m), rohelisega I boniteediklassi (H₁₀₀ = 29,5 m) ning punasega Ib boniteediklassi (H₁₀₀ = 37,5 m) kohta järgmised:

Enamikul juhtudel annab uus mudel suurema puistu tagavara juurdekasvu kui praegu kehtiv mudel;
Männi, kase, haava, sanglepa ja halli lepa puhul on eriti suur erinevus noores metsas, kus on suur täius (joonistel vanus 20 aastat ja täius 120%);
Männi puhul on viienda boniteedi puistu juurdekasv nooremates suurema täiusega metsas uue mudeli puhul väiksem kui praegu kehtiva mudeli järgi;
Kehvemate boniteetide puhul on mudelite erinevus väiksem kui viljakamatel metsamaadel;
Vanemates puistutes on suurem erinevus männi, kuuse ja haava puhul;
Vanemates kaasikutes ja lepikutes on puistu juurdekasvude vahe väike;
Noortes madala täiusega haavikutes ja hall-lepikutes on uuema mudeliga väiksem puistu juurdekasv, kui praegu kehtiva mudeliga arvutades;
Vanemates, suure täiusega kuusikutes on uue mudeliga arvutatud puistu juurdekasv palju suurem praegu kehtiva mudeliga arvutatust, eriti parema boniteediga puistutes.

Käesolevas töös loodud mudeliga ja „Metsa korraldamise juhendis“ (2009) oleva puistu tagavara juurdekasvu valemiga arvutatud puistu tagavara juurdekasvu erinevused puuliikide ja boniteetide kaupa: violetne – H₁₀₀ = 9,5 m, sinine – H₁₀₀ = 21,5 m, roheline – H₁₀₀ = 29,5 m, punane – H₁₀₀ = 37,5 m.
Differences in stand volume growth calculated with the model created in this study and the stand volume growth formula by the Forest Management Instruction (2009) by tree species and site index: violet – H₁₀₀ = 9.5 m, blue – H₁₀₀ = 21.5 m, green – H₁₀₀ = 29.5 m, red – H₁₀₀ = 37.5 m.

Mudelite rakendatavus ja andmete erisusest tulenevad probleemid

Siinses artiklis loodud mudel ei pruugi olla ideaalne kogu Eesti metsa jaoks, sest see põhineb Eesti metsa kasvukäigu püsiproovitükkide andmetel. Puistu kasvukäigu püsiproovitükkide puistute eelvalikul järgiti printsiipi, et proovitükkide arv igas metsatüübi rühmas oleks võrdeline selle metsatüübi pindalaga riigimetsas ja selle metsakasvukoha metsamaa hektarihinna korrutisega (Kiviste & Hordo, 2002). Rõhuv enamik proovitükkidest on rajatud riigimetsa, mis keskmiselt on enam hooldatud kui eramets, mistõttu andmestikus on suhteliselt suur osatähtsus majanduslikult eelistatumatel puuliikidel (mänd, kuusk, kask) ja ühevanuselistel puhtpuistutel.

Puistu rinnaspindala juurdekasv kasvukäigu püsiproovitükkidel on tõenäoliselt üsna täpne, sest viieaastase perioodi järel kordusmõõdeti proovitükkidel samade puude kahes suunas mõõdetud diameetreid. Andmestikus leidus siiski erindeid (suuri mõõtmiserinevusi), millest kontrollimisel selgus, et enamik neist ei olnud mõõtmisvead, vaid looduslikud iseärasused (Hordo et al., 2008), mis mõnel väiksemal proovitükil võis mõjutada tulemust. Loodud rinnaspindala juurdekasvu mudelit on võimalik kasutada puistu tagavara juurdekasvu arvutamiseks ka siis, kui muudetaks (parendataks) vormiarvu ja kõrguse kasvu mudeleid. Küll aga boniteedi ja/või täiuse normatiivi muutmise korral on vaja uuendada ka rinnaspindala juurdekasvu mudel.

Rinnaspindala juurdekasvu käsitlemine pideva funktsiooni tuletisena ja puistu kasvukäigu modelleerimine diferentsiaalvõrrandite süsteemina, mis pretendeerib ühtlasi bioloogilisele interpretatsioonile, on kõige enam välja arendatud Oscar García töödes (García, 2011; García et al., 2011). Siiski eelistatakse paljudes uurimustes (Allen II et al., 2020; Sun et al., 2007; Zhang et al., 2010) puistu tunnuste juurdekasvu modelleerimiseks üldistatud algebralise diferentsvõrrandi lähenemisviisi (Cieszewski & Bailey, 2000).

Täiuse arvutamise aluseks on standardtabelid („Metsa korraldamise juhend„, 2023: lisa 13), mis põhinevad normaalpuistu tagavara ja ristlõikepindala sõltuvusel ainult puistuelemendi keskmisest kõrgusest. Artur Nilson (2014) on kirjutanud, et täiuse arvutamine on tuginenud juba ca 100 aastat Eichhorni reeglile ning seda on metsameestele õpetatud ilma kriitikata ja alternatiivideta. Nilson (2014) toob peamiseks puuduseks boniteedi mittearvestamise, aga kritiseerib ka normaalpuistu tagavara mudelikuju. Alternatiivina pakub ta välja mudeleid, kus on kasutatud puistu tiheduse kirjeldamiseks piirhõreduse või ristlõikepindala seost puistu rinnasdiameetrist ja boniteedist.

Arvestades praegu kehtiva täiuse arvutamise puuduseid, on tulevikus otstarbekas ka üle vaadata metsanduse normatiivid, mis kirjeldavad puistu tihedust. Selles töös loodud mudeli puhul toob see kaasa parandi K uue modelleerimise ning põhimudelit selle puhul ei ole vaja uuesti luua.

Nii varasemas (Kohava, 1992) kui ka siinses töös esitatud mudelis arvutatakse segapuistu juurdekasv osapuuliikide juurdekasvude kaalutud keskmisena, millega eeldatakse juurdekasvus erinevate puuliikide vahelise olulise koosmõju puudumist. Soome ja Rootsi metsa kasvumudelite MOTTI ja HEUREKA testimine suure hulga segametsa katsealade andmetel näitasid, et nende mudelite vähemalt lühiajalised prognoosid osutusid headeks puuliikide erinevate koosseisude puhul (Aldea et al., 2023), mis kaudselt õigustab meie mudeli lihtsustust. Lisaks tehakse igal aastal Eesti metsa kasvukäigu proovitükkidel uusi mõõtmisi ning oleks otstarbekas suurema arvu mõõtmisandmete korral luua uus, korrigeeritud puistu juurdekasvu mudel.

Tänuavaldus

Eesti metsa kasvukäigu püsiproovitükkide võrgustiku rajamist ning kordusmõõtmisi on rahastanud SA Keskkonnainvesteeringute Keskus. Täname kõiki proovitükkide mõõtjaid. Artikkel valmis Riigimetsa Majandamise Keskuse ja Eesti Maaülikooli vahel sõlmitud lepingu „Puistuelemendi tagavara juurdekasvu mudeli väljatöötamine Eesti metsa kasvukäigu püsiproovitükkide võrgustiku and- metel ja selle kohta aruande esitamine“ tulemusena.

Summary

The stand volume and increment has been a topical issue in forest science and forest policy for over a hundred years. In today’s Estonia, representatives of various fields and ordinary people have increasingly started taking interest in this topic. It is important that people have a common understanding of different terms. In this work, the growth of the stand basal area was modelled as the sum of the net production of the basal area and basal area of dead trees (stand mortality) between the two consecutive measurements (Gadow et al., 2021; Pretzsch, 2009; Vaus, 2005).

The stand growth model must be usable in practice, i.e. traditional forest inventory data must been used as input (Peuhkurinen et al., 2008). The system of forest management models is constantly updated (Bravo et al., 2019; Kiviste et al., 2015). Therefore, it is justified to prefer the formula for the increase in the basal area of the stand over the increase in the volume of the stand (Kärki et al., 2000). The calculation of the increment of basal area is based on the measurement of the diameters but calculation of the stand volume growth includes additionally the height and stem taper.

The aim of this work was to create a new increment model of stand element for the most common tree species in Estonia based on the measurement data of individual trees. In this research, the network of permanent forest growth course sample plots (ENFRP) (Kiviste et al., 2015) from 1995 to 2020 has been used to create models according to the years of measurements. The data of the ENFRP network cover the growth site types and age classes of the stands of the majority of Estonian forests, while the locations of the plots in the stands have been chosen randomly. As a rule, the radius of ENFRP plots is 15, 20 or 25 m to ensure at least 100 first layer trees per sample plot.

Since repeated measurements of plots may not take place on the same date as the first measurement, the date of the measurement was also taken into account when calculating the age of the stand. Models were used to calculate the seasonality of growth for stand age, the development of which is described in a recently published article (Kiviste et al., 2022).

Since the aim was to model the gross increment of the basal area of the stand, it proved important to use the data of all trees that were alive at the beginning of the measurement period. It is not clear from the measurement data when the trees died. Thus, it was assumed that the mortality of trees during the study period was uniform in time, and trees were grown until the middle of the evaluated period. A methodology was used to calculate the diameter increment of dead trees, according to which the growth of thinner trees was lower and that of thicker trees it was larger. For this purpose, a methodology was created, in which a formula describing the growth of the diameter of all living trees of the first front is the first layer (equation 1). Then the relation of the parameters a and b of Formula 1 to the diameter is found (equation 5 and equation 6). Using the age of the tree that was halfway through the growing season in Formula 1 and parameters a and b, the last diameter of all the dead trees was found.

The basic formula (equation 9) for the calculation of the stand increment was based on the first year’s increase in the basal area of the first front (equation 7). When creating stand growth models, stands were selected with a height of at least 6 m and a basal area ratio of G₁ (basal area of first layer)/G_sum (total basal area of all layers of alive trees) of at least 0.95. Those pairs of measurements with a thinning volume between two measurements was greater than 5%, were excluded too. When building the Scots pine, Norway spruce and common alder models, only the sample plots were used, where the corresponding tree species was also the majority tree species. In the case of birch (silver and downy), European aspen and black alder, secondary tree species were also used in the modelling. The minimum proportion of those tree species was 10%. The data of the adjacent tree species were used for these tree species because there would have been too little data. There were too few data for common alder as well, but the increments of alder stand elements as the majority tree species and secondary tree species turned out to be too different. A total of 809 measurement pairs were used for Scots pine, 314 for Norway spruce, 373 for silver and downy birch, 273 for European aspen, 129 for black alder and 19 for common alder (Tables 1–6, respectively).

Before modelling, the measured increments of the basal areas of the stand elements of the first layer were divided by the proportion of the basal area of the stand element. The shape of the formula for the growth of the cross-sectional area of the stand was chosen such that the starting point would be (0.0), and as the age increases, the function would initially be increasing, later decreasing, and the function would not give a negative value in any case if the age were positive. A formula form corresponding to these characteristics and suitable for this work is equation 9, the properties of which are:

Starts at point 0.0;
The position of the maximum point relative to the x-axis (stand age) is described by the quotient of the estimates of parameters b₁ and b₂;
The maximum increment value depends on the quotient of the estimates of parameters b₀ and b₁;
The rate of decrease in the increment after the extremum point depends on the value of the estimate of the parameter b₁;
As age approaches infinity, the increment approaches 0.

Then the relationship between the parameters b₀, b₁ and b₂ of the main formula (equation 9) and the site index H₁₀₀ (equation 11, 15, 16, 18, 19, 20) was found. The indices of the parameters of the formula include tree species: “MA” – Scots pine; “KU” – Norway spruce; “KS” – downy and silver birch; “HB” – European aspen; “LM” – black alder and “LV” – common alder).

After estimating the parameters of the basic formula 9, the estimation of the correction K (actual basal area and the quotient of the increase in the basal area of the stand calculated by equation 9) was performed for each tree species. The need for correction K arises from the fact that the density of the stands is different, but the stand density is not used in the basic formula 9. Since the original data may not have the same proportion of stand density in all age classes and all site index classes, the need to include stand age and site index was checked when modelling the correction K. As a result, the shape of the correction K formula is somewhat different for each tree species (equation 21 to equation 24).

The visual results of the error analysis of the formula (measured stand growth minus model-calculated stand growth) are presented in the appendices: in Appendices 1–3 there is Scots pine, 4–6 Norway spruce, 7–9 silver and downy birch, 10–12 European aspen, 13–15 black alder and 16–18 common alder model errors compared with site index H₁₀₀, stand age and stand density, respectively.

Since the current stand growth model intended for use in Estonian conditions (Metsa korraldamise juhend, 2023: Appendix 12) is not compiled based on data from Estonian permanent sample plots, there was reason to believe that the results calculated with the model may give too large an error compared to reality. Figure 24 shows the differences in the calculation results of the stand volume growth by tree species (stand volume growth calculated with the model created in this work minus the currently valid stand volume growth): violet – H₁₀₀ = 9.5 m, blue – H₁₀₀ = 21.5 m, green – H₁₀₀ = 29.5 m, red – H₁₀₀ = 37.5 m.

The model created in this article is not ideal, but it is created based on the data of ENFRP. Since the model is primarily based on the growth of the basal area of the stand, the advantage of this model is that in the case of improving-changing the models of the stem taper or height growth, it is not necessary to calculate new parameters for the model. At the same time, the model uses site index H₁₀₀ and density, for which it is necessary to revise the increment model of the stand if the site index or density calculation models are changed.

The calculation of stand density is based on the dependence of the basal area of the normal stand only on the average height of the stand element. Artur Nilson (2014) has written that the calculation of fullness has been based on Eichhorn’s rule for about 100 years and it has been taught to foresters without criticism or alternatives. Nilson (2014) cites the failure to consider the site index as the main drawback, but also criticizes the model shape of the normal stand volume. As an alternative, he proposes models where the relationship between the basal area and the breast diameter and site index of the stand is used to describe the density of the stand. Considering the shortcomings of the currently valid stand density calculation, it will also be expedient to revise the forestry standards that describe the stand density in the future. In the case of the model created in this paper, this leads to a new modelling of the correction K, and it is not necessary to recreate the basic model (equation 9).

Both in the earlier model (Kohava, 1992) and in the model presented in this work, the growth of a mixed stand is calculated as a weighted average of the growth of individual tree species, with which it is assumed that there is no significant interaction between different tree species in the growth. The testing of Finnish and Swedish forest growth models MOTTI and HEUREKA on the data of a large number of mixed forest experimental areas showed that at least the short-term predictions of these models turned out to be good for different compositions of tree species (Aldea et al., 2023), which indirectly justifies the simplification of our model.

For users of the formula, the algorithm for the basal area increment of the stand, together with parameter estimates, is given in Appendix 19.

In addition, every year, new measurements are made on ENFRP sample plots, and it would be expedient to create a new, corrected stand growth model with a larger number of measurement data after some years.

Appendix

Appendix 1.

Prediction error of Scots pine stand basal area growth model depending on site index.

Männikute rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 2.

Prediction error of Scots pine stand basal area growth model depending on stand age.

Männikute rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 3.

Prediction error of Scots pine stand basal area growth model depending on stand relative density.

Männikute rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 4.

Prediction error of Norway spruce stand basal area growth model depending on site index.

Kuusikute rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 5.

Prediction error of Norway spruce stand basal area growth model depending on stand age.

Kuusikute rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 6.

Prediction error of Norway spruce stand basal area growth model depending on stand relative density.

Kuusikute rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 7.

Prediction error of birch stand basal area growth model depending on site index.

Kase rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 8.

Prediction error of birch stand basal area growth model depending on stand age.

Kase rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 9.

Prediction error of birch stand basal area growth model depending on stand relative density.

Kase rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 10.

Prediction error of aspen stand basal area growth model depending on site index.

Haava rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 11.

Prediction error of aspen stand basal area growth model depending on stand age.

Haava rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 12.

Prediction error of aspen stand basal area growth model depending on stand relative density.

Haava rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 13.

Prediction error of common alder stand basal area growth model depending on site index.

Sanglepa rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 14.

Prediction error of common alder stand basal area growth model depending on stand age.

Sanglepa rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 15.

Prediction error of common alder stand basal area growth model depending on stand relative density.

Sanglepa rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 16.

Prediction error of grey alder stand basal area growth model depending on site index.

Halli lepa rinnaspindala juurdekasvu vea sõltuvus boniteedist.

Appendix 17.

Prediction error of grey alder stand basal area growth model depending on stand age.

Halli lepa rinnaspindala juurdekasvu vea sõltuvus vanusest.

Appendix 18.

Prediction error of grey alder stand basal area growth model depending on stand relative density.

Halli lepa rinnaspindala juurdekasvu vea sõltuvus täiusest.

Appendix 19.

Stand basal area growth model algorithm.

Lisa 19. Mudeli kasutamine.

Valemi kasutamisel arvutatakse esmalt parameetrid b₀, b₁, b₂ ja parand K:

Mänd (pine):
$\begin{array}{l} b_{0, M A} & = - 0, 0012210 \cdot H_{100}^{2} + 0, 158686 \cdot H_{100} - 5, 9781; \\ b_{1, M A} & = - 0, 0003358 \cdot H_{100}^{2} - 0, 001375 \cdot H_{100} + 2, 1019; \\ b_{2, M A} & = \frac{b_{1, M A}}{0, 0410 \cdot H_{100} - 7, 081}; \\ K_{M A} & = {(H_{100})}^{- 0, 1216} \cdot {(\frac{A}{A + 1})}^{- 11, 5142} \cdot {(\frac{T}{T + 10})}^{0, 6164} \cdot {(\frac{G_{1}}{G_{s u m}})}^{0, 7270} . \end{array}$
Kuusk (spruce):
$\begin{array}{l} b_{0}_{, K U} & = 0, 0018364 \cdot H_{100}^{2} + 0, 010248 \cdot H_{100} - 4, 6873; \\ b_{1}_{, K U} & = - 0, 0022301 \cdot H_{100}^{2} + 0, 095526 \cdot H_{100} + 1, 0537; \\ b_{2}_{, K U} & = \frac{b_{1, K U}}{0, 0556 \cdot H_{100} - 7, 653}; \\ K_{K U} & = {(\frac{H_{100}}{H_{100} + 10})}^{- 3, 0160} \cdot {(\frac{A}{A + 1})}^{- 7, 3726} \cdot {(\frac{T}{T + 1000})}^{0, 4279} \cdot {(\frac{G_{1}}{G_{s u m}})}^{0, 4977} . \end{array}$
Kask (birch):
$\begin{array}{l} b_{0}_{, K S} & = 5, 2140 \cdot \ln (H_{100}) - 20, 3850; \\ b_{1}_{, K S} & = - 2, 2018 \cdot \ln (H_{100}) + 9, 7701; \\ b_{2}_{, K S} & = \frac{b_{1, K S}}{0, 05645 \cdot H_{100} - 6, 889}; \\ K_{K S} & = {(\frac{H_{100}}{H_{100} + 10})}^{- 0, 9095} \cdot {(\frac{A}{A + 1})}^{- 2, 7151} \cdot {(\frac{T}{T + 10})}^{2, 6006} \cdot {(\frac{G_{1}}{G_{s u m}})}^{1, 9882} . \end{array}$
Haab (aspen):
$\begin{array}{l} b_{0}_{, H B} & = 1, 9202 \cdot \ln (H_{100}) - 8, 4410; \\ b_{1}_{, H B} & = - 0, 4439 \cdot \ln (H_{100}) + 2, 9609; \\ b_{2}_{, H B} & = \frac{b_{1, H B}}{0, 02584 \cdot H_{100} - 6, 5293}; \\ K_{H B} & = {(\frac{H_{100}}{H_{100} + 10})}^{- 2, 1190} \cdot {(\frac{T}{T + 10})}^{5, 9814} . \end{array}$
Sanglepp (black alder):
$\begin{array}{l} b_{0}_{, L M} & = 1, 0000 \cdot \ln (H_{100}) - 4, 7467; \\ b_{1}_{, L M} & = 1, 2460; \\ b_{2}_{, L M} & = - 0, 2510; \\ K_{L M} & = {(\frac{H_{100}}{H_{100} + 10})}^{- 2, 6860} \cdot {(\frac{A}{A + 1})}^{- 4, 8688} \cdot {(\frac{T}{T + 1000})}^{0, 4772} \cdot {(\frac{G_{1}}{G_{s u m}})}^{0, 4716} . \end{array}$
Hall lepp (grey alder):
$\begin{array}{l} b_{0}_{, L V} & = 1, 0000 \cdot \ln (H_{100}) - 22, 7712; \\ b_{1}_{, L V} & = 11, 9181; \\ b_{2}_{, L V} & = - 1, 8336; \\ K_{L V} & = 0, 01158 \cdot T . \end{array}$

Parandi K arvutamisel on kaks erandit (seda ainult K arvutamisel ehk parameetrite b₀, b₁ ja b₂ arvutamisel alam-ega ülempiire ei ole):

parandi K arvutamisel on kuuse puhul alumiseks boniteedi (H₁₀₀) piiriks 9,5 ehk väiksema väärtusega H₁₀₀ puhul kasutatakse väärtust 9,5;
männi puhul on K arvutamisel alumiseks vanuse väärtuseks 10 aastat ehk kui männi element on sellest noorem, siis kasutatakse K arvutamisel vanuse väärtusena kümmet.

Pärast mudeli parameetrite ja parandi leidmist arvutatakse puistu ristlõikepindala juurdekasv järgmiselt:

z G = K \cdot e^{b_{0} + b_{1} \cdot \ln A + b_{2} \cdot {(\ln A)}^{2}} .

Mudeliga arvutatud puistu ristlõikepindala juurdekasvu kaudu puistu tagavara juurdekasvu arvutamiseks tuleb arvestada metsa kõrguse juurdekasvuga ja kasutada puistu vormiarvu mudelit. Puistuelemendi kõrguse kasvu arvestamiseks on võimalik kasutajal valida kahe Eesti mudeli vahel. Esimene variant on kasutada „Metsa korraldamise juhendi“ (2023) lisas 10 toodud boniteedi H₁₀₀ arvutamisvalemit puistuelemendi kõrguse leidmiseks aasta pärast järgmiselt:

\begin{array}{l} H_{50} = \frac{H_{100} \cdot {1 + α \cdot [{(\frac{50}{100})}^{c} - 1]}}{{1 - β \cdot H_{100} \cdot [{(\frac{50}{100})}^{c} - 1]}}; \\ H_{A + 1} = \frac{H_{50}}{{1 + (α + β \cdot H_{50}) \cdot [{(\frac{50}{A + 1})}^{c} - 1]}}, \end{array}

kus H₁₀₀ on boniteet (m), A on puistu vanus, a, β ja c on valemi konstandid.

Teine võimalus on kasutada Kiviste (Kiviste & Kiviste, 2009) moodpuistute diferentsmudelit:

\begin{array}{l} β H = c_{4} - 493 \cdot \ln (O H O R + 1) + 1355 \cdot k \cdot k P; \\ d H = β H / 50^{c_{1}}; \\ r H = \sqrt{{(H - d H)}^{2} + 4 \cdot β H \cdot H / A^{c_{1}}}; \\ H_{A + 1} = (H + d H + r H) / (2 + 4 \cdot β H \cdot {(A + 1)}^{(- c_{1})} / (H - d H + r H)), \end{array}

kus H_A+1 on puistuelemendi kõrgus aasta pärast, m; βH, dH, rH – abimuutujad; c₁, c₄ – puuliigipõhised konstandid; kP – okaspuude puhul 1, lehtpuude korral 0; OHOR – kõduhorisondi tüsedus, cm; k – kui kultuur, siis 1, muul juhul 0; A – puistuelemendi vanus, a; H – puistuelemendi kõrgus, m.

Pärast kõrguse arvutamist on vaja arvutada kaks vormiarvu, üks vanuses A kõrgusega H ning teine vanuses A+1 kõrgusega H_A+1. Vormiarvu (F) arvutamiseks on võimalik kasutada kirjandusallikates leitavatest sobilikke mudeleid või „Metsa korraldamise juhendi“ (2023) lisas 11 toodud mudelit.

Puistu tagavara juurdekasvu saab seega arvutada puistu ristlõikepindala juurdekasvu kaudu valemiga:

z M = F_{A + 1} \cdot (G + K F_{G} \cdot z G) \cdot H_{A + 1} - F \cdot G \cdot H,

kus zM on puistuelemendi tagavara juurdekasv, m³/ha/aasta; F_A+1 – puistuelemendi vormiarv aasta pärast; G – puistuelemendi ristlõikepindala enne kasvatamist, m²/ha; KF_G – puistuelemendi G osakaal (G_element/G_{I rinne}); zG – puistuelemendi G juurdekasv juhul, kui puistuelemendi koosseis oleks 100%, m²/ha/aasta; F – puistuelemendi vormiarv enne kasvatamist; H – puistuelemendi kõrgus enne kasvatamist, m.

The model of stand basal area gross growth on the data of the Estonian Network of Forest Research Plots

Full Article

Sissejuhatus

Andmed ja metoodika

1

2

3

4

5

6

7

8

9

10

Tabel 1.

Tabel 2.

Tabel 3.

Tabel 4.

Tabel 5.

Tabel 6.

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

Tulemused ja arutelu

Puistu ristlõikepindala juurdekasvu mudeli parameetrid ja headus

Tabel 7.

Tabel 8.

Tabel 9.

Tabel 10.

Joonis 1.

Joonis 2.

Tabel 11.

Joonis 3.

Joonis 4.

Joonis 5.

Joonis 6.

Joonis 7.

Joonis 8.

Joonis 9.

Joonis 10.

Joonis 11.

Joonis 12.

Joonis 13.

Joonis 14.

Joonis 15.

Joonis 16.

Joonis 17.

Joonis 18.

Joonis 19.

Joonis 20.

Joonis 21.

Joonis 22.

Joonis 23.