215 t u t k i m u s a r t i k k e l i Metsätieteen aikakauskirja Annika Kangas Jouni Siipilehto Jouni Siipilehto ja Annika Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta suomalaisissa talousmetsissä Siipilehto, J. & Kangas, A. 2015. Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta suomalaisissa talousmetsissä. Metsätieteen aikakauskirja 4/2015: 215–236. Tutkimuksessa tarkasteltiin Näslundin pituuskäyrän matemaattisia ominaisuuksia ja laadittiin malliperhe Näslundin pituuskäyrän ennustamiseksi. Aineistona käytettiin taimikoiden TINKA- ja varttuneiden puustojen INKA-kokeita. Pituuskäyrät sovitettiin männylle, kuuselle ja lehtipuille Näslundin yhtälön linearisoidussa muodossa. Ennustemallit laadittiin havumetsissä pääpuulajil- le, mutta lehtipuiden aineisto koostui sekä koivikoista, että sekapuuna kasvaneesta koivusta ja vähäisestä määrästä muuta lehtipuustoa. Näslundin pituuskäyrää ennustettiin sekä nuorten että varttuneiden metsiköiden tyypillisillä puustotunnusten yhdistelmillä ja lisäksi valtapituudella ja puuston tilavuustunnuksilla. Ennustemallit laadittiin lineaarisina sekamalleina, jotta mallit ovat lokalisoitavissa lineaarisen ennustamisen teorian mukaisesti estimoimalla satunnaisosat puutason mittauksilla. Jäännöshajonta oli heteroskedastinen metsikön kehitysvaiheen, keskimääräisen so- lakkuuden sekä lämpösumman suhteen metsiköiden välillä ja se huomioitiin varianssifunktion avulla. Jäännösvirheen hajontaa voitiin hyödyntää pituuden satunnaisvaihtelun luomiseksi. Pääosa laadituista malleista on tarkoitettu vaihtoehdoksi tilanteessa, jossa puutason tietoja ei ole. Mallin kiinteän osan avulla saatiin varsin luotettava kuva puiden läpimitan ja pituuden riippuvuudesta. Laaditut mallit on helppo ottaa käyttöön esim. erilaisissa metsätalouden suunnittelujärjestelmissä. Avainsanat: ennustaminen, kalibrointi, puustotunnukset, sekamalli Yhteystiedot: Luonnonvarakeskus (Luke), Vantaa & Joensuu Sähköposti jouni.siipilehto@luke.fi Hyväksytty 15.12.2015 Saatavana http://www.luke.fi/aikakauskirja/full/ff15/ff154215.pdf 216 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli 1 Johdanto Puun pituuden ja läpimitan välinen riippuvuus on kautta aikain pysynyt metsien mallintami- sessa mielenkiinnon kohteena. On olemassa joukko tutkimuksia, joissa useita vaihtoehtoisia funktioita verrataan keskenään läpimitan ja pituuden riippu- vuuden kuvaamiseksi. Esimerkkinä tällaisista mai- nittakoon Prodan (1965), Curtis (1967), Arabatzis ja Burkhart (1992), Elfving ja Kiviste (1997), Zhang (1997), Fang ja Bailey (1998), Newton ja Ampon- sah (2007), Leduc ja Goelz (2009) ja Gómez-García ym. (2014) ja testattuja funktioita oli esimerkiksi Korf, Gompertz, Schumacher, Chapman-Richards ja Weibull. Minkään tietyn funktion ei ole todettu olevan selvästi muita parempi (ks. Mehtätalo 2004, Mehtätalo ym. 2015). Tyypillisesti Näslundin (1936) pituuskäyrää ei näissä vertailuissa ollut mukana tai se oli mukana, mutta sen alkuperää ei tunnistettu (esim. Fang ja Bailey 1998, Staudhammer ja Le- May 2000, Leduc ja Goelz 2009). Pituuskäyräksi suositellaan sigmoidia funktiota, jolla on sekä kään- nepiste että tunnettu asymptootti (esim. Yuancai ja Paresol 2001). Näslundin pituuskäyrän käännepiste ja asymptootti tunnetaan, mutta asymptootti ei ole suoraan yksittäinen parametri, kuten muissa edellä mainituissa funktiossa. Aivan tuoreessa tutkimuksessa Näslundin pituus- käyrä oli mukana vaihtoehtoisten mallien vertaillus- sa ja se pärjäsi erinomaisesti (Mehtätalo ym. 2015). Näslundin pituuskäyrä on varsin helppokäyttöinen muun muassa siksi, että siinä on vain kaksi esti- moitavaa parametria. Tyypillisesti Näslundin pi- tuuskäyrää käytetään koepuiden läpimitta- ja pi- tuushavaintojen tasoittamiseksi, jotta koepuutiedot voidaan yleistää lukupuille tasoituskäyrältä. Metsän- tutkimuslaitoksella kehitetyssä KPL-laskentaohjel- massa Näslundin pituuskäyrä on oletusarvona ja sen vaihtoehtona on Prodanin (1965) pituuskäyrä (ks. Heinonen 1994). Myös tietyt simulointiohjelmat, kuten MOTTI Suomessa ja HEUREKA Ruotsissa soveltavat Näslundin pituuskäyrää yhtenä vaihtoeh- tona (Fahlvik ym. 2014, Siipilehto ym. 2014). Pituuskäyrän ennustemalleja tarvitaan tyypillisesti silloin, kun puuston rakennetta halutaan ennustaa puutasolla. Puiden läpimitat poimitaan ennuste- tusta läpimittajakaumasta ja vastaava puun pituus ennustetun pituuskäyrän avulla. Siipilehdon (1999) sekametsäaineistoista (ks. Mielikäinen 1980, 1985) laatimia Näslundin pituuskäyrän ennustemalleja on käytetty Suomessa aika yleisesti (esim. Maltamo ym. 2006, Kuusisto ja Kangas 2008, Mustonen ym. 2008, Kangas ym. 2010, 2012, Muinonen ym. 2013). Siipilehto (2011b) havaitsi niiden tuottavan harhaisia tuloksia mallin laadinta-aineiston ulkopuo- lella, erityisesti nuorissa metsissä. Tämä systemaat- tinen virhe johtui pääasiassa siitä, että pituuskäyrän parametrien riippuvuus oli kuvattu suoraan keski- läpimitan ja keskipituuden lineaarisena funktiona (Siipilehto 1999). Myöhemmin Siipilehto (2011a) linearisoi parametrien ja puustotunnusten välisen riippuvuuden ottamalla logaritmit molemmista. Lineaarisen ennustamisen menetelmällä Näslun- din pituuskäyrän parametrit voidaan lokalisoida (kalibroida) millä tahansa tunnetuilla puustotun- nuksilla (Siipilehto 2011a) tai puutason mallissa perustuen koepuuotantaan (Kangas ja Maltamo 2002, Kinnunen ym. 2007, Vastaranta ym. 2010, Schmidt ym. 2011). Menetelmän soveltaminen vaa- tii matriisilaskentaa parhaan lineaarisen ennustimen (best linear unbiased predictor) eli BLUP-estimaatin laskemiseksi (ks. Lappi 1993, Kangas ja Maltamo 2002). Siitä syystä niiden käyttöönotto ei ole aivan yksinkertaisen suoraviivaista vaan vaatii tavallisia regressiomalleja enemmän esimerkiksi ohjelmoin- titaitoja. Hyvin yleisen pituusmallin lokalisointi yhdelläkin mittaushavainnolla korjaa ennustetta voimakkaasti havainnon suuntaan (ks. Lappi 1991, Mehtätalo 2005). Siten täysin satunnaiseen mittaus- havaintoon kohdistuu pieni riski. Jos satunnainen puu osuu läpimitan suhteen erityisen lyhyeen tai pitkään puuhun, lokalisoitu malli automaattisesti aliarvioi tai jälkimmäisessä tapauksessa yliarvioi metsikön pituuskäyrää. Useampaa mittaushavain- toa käytettäessä riski pienenee. Puustotunnuksilla lokalisoitava Näslundin pituuskäyrän parametrien ennustemalli osoittautui oivaksi vaihtoehdoksi Korf- funktion (Mehtätalo 2005) tai Vetlheimin (1987) pituusmallin rinnalle (ks. Siipilehto 2011b, s. 29). Koska Siipilehto (2011b) ennusti pituuden sijaan pituuskäyrän parametreja, eivät kyseisen mallin virhetermit kuvanneet mitään sen luotettavuudesta itse puun pituuden ennusteena. Tässä raportissa esitellään Näslundin pituuskäy- rän yleisiä ominaisuuksia ja regressiomallinnuksen 217 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… vaihtoehtoja. Vaihtoehtoisina menetelminä tarkas- tellaan regressiomalleja suoraan puun pituudelle sekä linearisoidulle muunnokselle että etukäteen sovitetuille Näslundin pituuskäyrän parametreille. Puukohtaiset mallit estimoidaan puun pituudelle epälineaarisena ja muunnokselle lineaarisena seka- mallina. Lopuksi laaditaan uusi malliperhe männyn, kuusen ja koivun pituuden ennustamiseksi. Tavoit- teena on mallien harhaton käyttäytyminen varttu- neista taimikoista aina uudistuskypsiin metsiköihin. Regressiomalleista esitetään useita vaihtoehtoja, joista voidaan valita puustotunnuksiltaan käyttä- jän aineistoihin ja suunnittelujärjestelmiin sopivin vaihtoehto. Vaihtoehtoisia puustotunnuksia tarvitaan esim. puuston kehitysvaiheen ja puuston käsittelyn mukaan. Niinpä mallien laadinnassa huomioidaan erot tyypillisissä puustotunnuksissa taimikoissa ja nuorissa metsissä verrattuna varttuneisiin metsiin tai leimikoiden tietoihin. 2 Aineisto ja menetelmät 2.1 Aineisto Aineisto edusti suomalaisia talousmetsiä. Se koostui valtakunnan metsien inventointien VMI7- ja VMI8- koealaverkkoon perustetuista varttuneempien puus- tojen INKA- ja taimikoiden TINKA-kokeista (Gus- tavsen ym. 1988). Kokeet on perustettu 1976–1986 välisenä aikana. INKA on mitattu neljään ja TINKA kolmeen kertaan. Metsikkökoeala koostuu kolmen ympyräkoealan rypäästä. Ympyräkoealojen koko valittiin puuston tiheyden mukaan siten, että kol- melta koealalta mitattiin yhteensä noin 120 puuta eteläisessä Suomessa ja noin 100 puuta Lapissa. INKA-kokeiden ympyräkoealan sisältä mitattiin pienempisäteinen koeala koepuita varten siten, että koepuukoeala edusti kolmannesta metsikkökoealas- ta. Koepuista mitattiin läpimitan lisäksi pituus ja esim. latvustunnuksia. TINKA-kokeiden kahdella ensimmäisellä mittauskerralla kaikista rinnankor- keuden ylittäneistä puista mitattiin sekä läpimitta että pituus. Koealat yhdistettiin edustamaan koko metsikköä, jolloin pituuskäyrän sovitukseen jäi kes- kimäärin 40 puuta eteläisen Suomen ja 33 puuta Lapin INKA-aineistossa. Koska koivumetsiköitä oli vähän ja koivu sekä muut lehtipuut esiintyvät pääasiassa sekapuustona havumetsissä, koostettiin lehtipuiden mallitusaineisto yhdistämällä koivikot (27 metsikköä) ja koivu- ja lehtisekapuustot (39 met- sikköä). Sekapuustoissa lehtipuita piti olla vähintään 10 kpl metsikössä. Lehtipuuryhmä koostui lähes yk- sinomaan hies- ja rauduskoivusta. Todellinen ikä oli aina määritetty metsikön pääpuulajista ja se rajoi- tettiin lehtisekapuustoissa 120 vuoteen. Metsikön todellinen ikä perustuu valtapuiden keskimääräiseen rinnankorkeusikään ja INKA-kokeiden kairauksista laadittuun talousmetsien ikälisäysmalliin (Siipilehto ja Huttunen 2015). Näin määritelty todellinen ikä oli pienempi kuin aiemmin yleisesti käytetty VMI- ikälisäykseen (Kuusela ja Salminen 1969) perus- tunut todellinen ikä. Alikasvos ja ylispuut jätettiin pituuskäyrän sovituksen ulkopuolelle. Taimikoiden TINKA-aineisto rajattiin keskipituudeltaan yli 4 m metsiköihin, koska käytännöllisesti katsoen kaikki puut ovat saavuttaneet tässä vaiheessa rinnankorke- uden (ks. Siipilehto 2011a) ja siten niille saadaan relevantit rinnankorkeuteen perustuvat puustotun- nukset. Lopullinen aineisto koostui 568 männiköstä, 214 kuusikosta sekä 66 koivikosta tai lehtisekapuus- tosta. Puukohtaisia havaintoja oli männyllä 32 037 kuusella 8 339 ja lehtipuista 1 676. 2.2 Näslundin pituuskäyrän yleisiä ominaisuuksia Näslund (1936) esitteli pituuskäyrän männylle käyt- täen toisen asteen yhtälöä. Kuitenkin kolmannen asteen yhtälö oli joustavampi kuvaamaan kuusen tyypillisesti laajaa läpimitta ja pituusvaihtelua (Pet- terson 1955, Vestjordet 1972, Siipilehto 1999). Siten Näslundin pituuskäyrän yleisempi muoto voidaan kirjoittaa yhtälön 1 mukaisesti: h = dm b0 + b1d( )m +1,3 (1) jossa h on puun pituus, d on puun läpimitta rinnan- korkeudelta (1,3 m), potenssi m = 2 valopuulajeille (esim. mänty ja koivu) ja m = 3 varjoa sietäville la- jeille (esim. kuusi), b0 ja b1 ovat estimoitavia pa- rametreja. Pituusmallissa on siten kaksi estimoi- tavaa parametria, koska potenssia m ei estimoida, vaan se valitaan etukäteen puulajin ominaisuuksiin 218 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli perustuen. Näslundin pituuskäyrän asymptootti saadaan b1 parametrista ja käyrän potenssista m ja käännepiste saadaan molempien parametrien avulla. Käännepis- te (d’’, h’’) saadaan toisen derivaatan nollakohdan positiivisesta juuresta. Asymptootti on (1 / b1 m + 1,3) ja käännepiste on: d’’ = (m – 1)b0 / 2b1 (2) ja h’’ = (b0(m + 1))–m(b0(m – 1) / b1)m + 1,3 (3) Kaavat 2 ja 3 sievenevät, kun sijoitetaan m:n arvoiksi 2 ja 3. Kun m = 2, d’’ = b0 / 2b1 ja h’’ = 1 / 9b1 2 (ks. Näslund 1936, s. 43), kun m = 3, d’’ = 2b0 / 2b1 ja h’’ = 1 / 8b1 3. Pituuskäyrän parametrit voidaan estimoida epäline- aarisella regressiolla suoraan yhtälöstä (1). Yleisem- min pituuskäyrän parametrit estimoidaan parametri- en suhteen linearisoidun yhtälön (4) avulla: y = d h −1,3( )m −1 = b0 + b1d + ε y (4) Koska pituudelle tehdään epälineaarinen muunnos, tulee se ottaa huomioon harhankorjauksena takaisin- muunnoksen yhteydessä. Taylorin sarjaan perustuva harhan korjaus on muotoa E(h) = g(y) + (1/2)g’’(y)s2, missä g on käänteisfunktio (ks. Lappi 1993, s. 91). Eerikäinen ja Korhonen (2001, s. 75) esittivät har- hakorjatun pituusennusteen Näslundin toisen asteen yhtälölle:  h = d 2 b0 + b1d( )2 + 3d 2 b0 + b1d( )4 σ 2 +1,3 (5) Vastaava harhakorjattu pituusennuste kuuselle Näslundin kolmannen asteen yhtälöstä saadaan seuraavasti. Ensin ratkaistaan käänteisfunktio, joka saadaan: y = d h −1,3( ) 1 3 ⇔ y h −1,3( ) 1 3 = d⇔ h −1,3( ) 1 3 = d y ⇔ h −1,3= d 3 y3 ⇔ h = d 3y−3 +1,3 (6) Käänteisfunktion ensimmäinen ja toinen derivaatta ovat: ′g y( ) = d 3 −3y( )−4 ja ′′g y( ) = d 3 −3( ) −4y( )−5 = 12d 3y−5 (7) Toisen derivaatan mukainen harhan korjaustermi kuuselle on siten: 1 2 12d 3 y5 σ 2 = 6d 3 b0 + b1d( )5 σ 2 (8) Täten kuusen lopullinen harhakorjattu pituusennuste saadaan:  h = d 3 b0 + b1d( )3 + 6d 3 b0 + b1d( )5 σ 2 +1,3 (9) Linearisoidun mallin etuna on se, että jäännösvirhe voidaan olettaa homogeeniseksi ja normaalijakau- tuneeksi eli yksittäisen metsikön sisällä varianssi on vakio (Näslund 1936, s. 52). Kun lineaarisen yhtälön (4) normaalijakautuneen virheen keskihajonta (sy) palautetaan alkuperäiseen skaalaan pituuden virhe- vaihteluksi (sh), on hajonta läpimitan ja pituuden funktio yhtälön (10) mukaisesti (Näslund 1936, s. 56). Yhtälö johdetaan Taylorin sarjan avulla (ks. Siipilehto 2000). Yhtälöä (10) tarvitaan, jos pituu- den odotusarvon lisäksi halutaan kuvata realistista metsikkökohtaista pituuden satunnaisvaihtelua li- nearisoitua yhtälöä (4) käytettäessä. sh = sy m ĥ −1,3( ) m+1 m ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ d (10) Jos pituuskäyrältä tunnetaan jokin piste, voidaan Näslundin pituuskäyrän toinen parametri ratkais- ta siten, että yhtälö kulkee kyseisen pisteen kautta (Siipilehto 1999). Parametreista b1 on osoittautu- nut parametria b0 voimakkaammin puustotunnus- ten kanssa korreloituneeksi. Oletetaan siis, että b1 ennustetaan puustotunnusten, kuten metsikön pohjapinta-alan mediaaniäpimitan (DGM) ja pi- tuuden (HGM) avulla. Jotta pituuskäyrä kulkisi pis- teen (DGM, HGM) kautta, ratkaistaan b0 yhtälöllä: b0 = DGM / (HGM – 1,3)(1/m) – b1 DGM. Teoreettisesti parempi menetelmä on laatia seka- malli ja lokalisoida pituuskäyrä koepuumittauksen tai -mittausten avulla (ks. Kangas ja Maltamo 2002, Eerikäinen 2009). Menetelmässä korjataan mallin kiinteän osan ennustetta estimoimalla mallin sa- tunnaisparametrit mittaustiedon avulla. Esimer- 219 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… kiksi Kangas ja Maltamo (2002) estimoivat yleisen Näslundin pituuskäyrän männylle koko aineiston yli siten, että malli sisälsi pelkästään kiinteät parametrit b0 ja b1 sekä satunnaisen koealavaikutuksen (uij), satunnaisen vakion (b0i) että satunnaisen kertoimen (b1i) metsikölle i ja koealalle j. Malli lokalisoitiin mediaanipuun läpimitan ja pituuden (DGM, HGM) avulla. Koska tällainen malli lokalisoituu varsin te- hokkaasti, laadittiin tässä tutkimuksessa vastaava malli männylle, kuuselle ja lehtipuille. Kuvassa 1 on esitetty Näslundin pituuskäyrän ominaisuuksia, kuten tyypillinen pituuskäyrän muoto männyllä ja varjoa sietävällä kuusella. Pi- tuuskäyrät vastasivat aineiston keskiarvoja (tau- lukko 1). Kuusella erottuu käyrän sigmoidi muoto, koska käyrän käännepiste (d’’ = 5,0 ja h’’ = 4,8) on kauempana kuin männyllä (d’’ = 0,14 ja h’’ = 1,31). Männyllä pituuskäyrän käännepiste oli tyypillisesti niin lähellä käyrän lähtöpistettä eli 1,3 m:n pituutta, että käytännössä voidaan puhua konkaavista män- nyn pituuskäyrästä. Männyn pituuskäyrä taipui voi- makkaammin kuin kuusen ja siten sen asymptootti (18,4) oli pienempi kuin kuusella (29,2). Tällaiset erot näyttävät olevan tyypillisiä valopuulajien ja varjoa sietävien lajien välillä (Temesgen ja Gadow 2004). 2.3 Näslundin pituuskäyrän sovittaminen ja ennustaminen Ensiksi pituuskäyrät sovitettiin vaihtoehtoisilla me- netelmillä männylle siten, että kussakin mallissa oli samat selittäjät ln(DGM) ja ln(HGM). Vaihtoehtoi- na tarkasteltiin mallien kiinteää osaa menetelmällä A) etukäteen sovitettujen pituuskäyrän parametrien ennustaminen B) lineaarisen muunnoksen (kaava 4) ennustaminen puutason sekamallilla ja C) suora pituuden ennustaminen epälineaarisesti (kaava 1), jossa mallinnettiin logaritmisia parametreja Mehtä- talon (2015) lmfor paketin HDnaslund4 funktiolla eli h = (d / exp(b0) + exp(b1)d)2 + 1,3. Näitä laadittuja mallivaihtoehtoja verrattiin myös Siipilehdon (1999) pituusmalliin. Siipilehto (1999) laati mallin suoraan b1 parametrille DGM ja HGM selittäjien avulla ja parametri b0 ratkaistiin siten, että käyrä kulki pis- teen (DGM, HGM) kautta. Sekä menetelmä A) että C) tuottavat ennusteet logaritmisille parametreille ja siksi ne saavat aina varmasti positiivisen arvon. Vaihtoehtojen tarkastelun jälkeen päädyttiin linea- risoidun yhtälön (4) käyttöön. Lopullista mallinnusta varten sovitus tehtiin line- aarisen yhtälön (4) mukaisella regressiomallilla erik- seen jokaiselle metsikölle ja mittauskerralle puulaji- ryhmittäin. Nämä alustavat sovitukset palvelivat pa- Kuva 1. Männyn ja kuusen aineiston keskiarvoja vastaavat pituus- käyrät, niiden asymptootit ja käännepisteet. Pituuskäyriin on lisätty ±2 kertaa pituuden keskihajonta (sh) yhtälön 10 mukaan. 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Pi tu us , m Läpimitta, cm Männyn ja kuusen pituuskäyrän ominaisuuksia Mänty Männyn asymptootti Männyn käännepiste Kuusi Kuusen asymptootti Kuusen käännepiste 220 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli rametrien sekä jäännöshajonnan ja puustotunnusten välisten riippuvuuksien tarkastelua ja linearisointia. Lopullisen mallin sovituksessa jäännösvirheen kes- kihajonta kuvattiin varianssifunktion avulla. Kes- kimääräiset pituuskäyrän parametrit, linearisoidun mallin jäännöshajonta ja niitä selittävät puustotun- nukset puulajeittain on kuvattuna taulukossa 1 vaih- teluväleineen. Metsikköön sovitetuille pituuskäyrän parametreil- le laadittiin alustavia regressiomalleja metsikön mit- tauskertakohtaisten puustotunnusten funktiona. Al- kuperäinen riippuvuus parametrien ja puustotunnus- ten, etenkin parametrin b1 ja dimensioiden D, DGM, DDOM, H, HGM ja HDOM välillä linearisoitui, kun molemmille tehtiin logaritmimuunnos (ks. Siipilehto 2011b, kuva 2). Alustavista parametrien ennustemal- leista saatiin lineaarisen sekamallin potentiaalinen rakenne (ja tarvittaessa alkuarvot epälineaarisen sekamallin sovitukseen). Sekamallin ensimmäi- sessä versiossa oli sekä metsikkö että mittauskerta satunnaisena vakiona ja satunnaisena kertoimena. Mittauskerran kohdalla näiden välinen korrelaatio oli mallin estimoinnin kannalta liian voimakas (≈1). Ensin mittauskerran vaikutus poistettiin satunnai- sesta kertoimesta. Seuraavan vaiheen sovituksessa satunnainen mittauskertakohtainen vakio osoittautui tarpeettomaksi (lähes nolla) ja mallin hyvyyttä mit- taava BIC (Bayesian information criterion) parani, kun sekin jätettiin mallista pois. Lopulliset lineaa- riset sekamallit olivat yleisessä muodossa: yij = xi0’b0 + b0i + (xi1’b1 + b1i)dij + eij (11) jossa y on yhtälön (4) mukainen muunnos, xi0 on vaihtoehtoisten mallien kiinteän osan selittäjät (esim. mittauskerta- ja metsikkökohtaiset N, D ja H) parametrille b0 ja b0 niitä vastaavat kertoimet, xi1 on vaihtoehtoiset selittäjät parametrille b1 ja b1 on niitä vastaavat kertoimet, b0i on satunnainen met- sikkötason i vakio ja b1i satunnainen metsikkötason kerroin ja dij on puun j läpimitta metsikössä i ja eij on puutason jäännösvirhe. Satunnaisvaikutukset ole- tetaan normaalijakautuneiksi eli (b0i, b1i)~N(0, D), jossa D sisältää satunnaisvaikutusten varianssit ja kovarianssin. Jos jäännösvirhe on heteroskedastinen, se korja- taan varianssifunktion avulla. Usein käytetty vaih- toehto on ns. potenssifunktio, jossa vapaasti esti- Taulukko 1. Keskimääräiset pituuskäyrän parametrit (b0, b1), linearisoidun mallin jäännöshajonta (sy) ja niitä selittävät puustotunnukset puulajeittain ja metsiköiden lukumäärä (n). Koivu/lehtipuu aineisto sisälsi 27 koivikkoa ja koivu- tai lehtisekapuustoa 39 metsiköstä. (Puustotunnukset: DDY = lämpösumma, T = todellinen ikä, G = metsikön pohjapinta-ala, N = runkoluku, D = keskiläpimitta, H = keskipituus, DGM = pohjapinta-alan mediaaniläpimitta, HGM = pohjapinta-alan mediaanipituus, DDOM = valtaläpimitta, HDOM = valtapituus, Vtot = kokonaisrunkotilavuus. Mänty n = 568 Kuusi n = 214 Koivu n = 66 Keskim. Min. Maks. Keskim. Min. Maks. Keskim. Min. Maks b0 1,195 0,432 2,499 1,635 0,698 2,492 0,898 0,434 1,767 b1 0,242 0,150 0,432 0,330 0,265 0,457 0,242 0,150 0,449 sy 0,238 0,081 0,889 0,227 0,117 0,622 0,169 0,063 0,461 DDY 984 674 1348 1110 746 1356 1032 710 1305 T 62,0 11,0 195,0 79,4 18,0 173,0 64,8 10,0 120,0 G 13,8 0,9 52,2 21,0 2,4 62,3 9,2 0,4 25,7 N 1153 118 4552 1088 203 3183 975 115 7294 D 13,0 4,1 34,4 15,3 4,7 34,1 10,8 3,3 29,5 H 10,6 4,0 27,3 13,0 4,1 27,9 11,3 3,5 26,5 DGM 16,1 5,0 37,0 19,7 6,2 36,8 13,4 3,9 32,0 HGM 12,2 4,4 28,1 16,1 5,2 29,0 12,9 4,1 28,3 DDOM 20,9 8,2 41,4 26,5 9,0 43,4 17,2 5,9 34,8 HDOM 13,5 4,8 28,6 18,7 6,4 30,2 14,1 4,6 29,3 Vtot 90,5 5,4 567,1 171,2 7,6 651,3 63,1 0,7 252,3 Tukki 40,1 0 478,3 88,6 0 480,9 19,4 0 203,4 Kuitu 44,2 1,5 180,7 76,8 2,1 309,3 36,8 0 176,1 221 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… moitava potenssi (p) voi kuvata selittävän muuttujan suhteen pienenevää varianssia (p < 0) tai kasvavaa varianssia (p > 0) (ks. Mehtätalo ym. 2015). Tässä tutkimuksessa yksittäisen metsikön sisällä varianssi oletettiin vakioksi, mutta eri metsiköiden välillä va- rianssit poikkesivat toisistaan. Siksi varianssifunktio kuvattiin lämpösumman ja puuston keskimääräisen solakkuuden vaikutuksena metsiköiden välisen vari- anssin vaihtelussa. Varianssi kasvoi lämpösumman pienetessä ja solakkuuden (H / D tai HGM / DGM) pienetessä. Solakkuus pienenee sekä puuston vart- tuessa että etelästä pohjoiseen, joten voidaan sanoa, että solakkuudessa on mukana kilpailun lisäksi metsikön kehitysvaihe ja maantieteellinen sijainti. Varianssifunktiota käytettäessä perinteisestä jään- nöshajontatermistä s(eij) tulee varianssifunktion skaalausparametri. Tässä tutkimuksessa varians- sifunktion selittäjät skaalattiin siten, että keskiar- voksi tuli likimain yksi. Tällaisia selittäjiä olivat ln(H / D + 2), ln(HGM / DGM + 2) ja (1000 / DDY). Siten mallin estimoitu varianssifunktion skaalaus- parametri kuvasi likimain aineiston keskimääräistä puutason jäännöshajontaa ja varianssifunktio korjasi sitä molempiin suuntiin. Varianssifunktio huomioiden lopullinen y:n jään- nöksen metsikkökohtainen keskihajonta saatiin so- lakkuudesta sy(Hk, Dk) = syln(Hk / Dk + 2)p, jossa sy on estimoitu varianssifunktion skaalausparametri ja p on varianssifunktion estimoitu potenssi ja in- deksi k viittaa joko taimikon tai varttuneen puuston keskitunnuksiin. Lämpösummasta laskettaessa se oli sy(DDY) = sy (1000 / DDY)p. Näiden termien avulla ennustettiin metsikkökohtainen keskihajonta muun- noksen y jäännösvirheelle, joka kaavaa 10 käyttäen voitiin muuttaa puun pituuden keskihajonnaksi (ks. kuva 3). Metsikkökohtaista keskihajontaa tarvittiin myös harhan korjauksen yhteydessä (kaavat 5 ja 9). Metsiköistä arvioitavat puustotunnukset muuttu- vat puuston kehitysvaiheen mukaan. Tiheystunnus on tyypillisesti runkoluku (N) nuorissa metsissä ja pohjapinta-ala (G) varttuneissa metsissä, kuten esim. SOLMU- ja VMI-kuviotiedoissa (ks. Sol- mun… 1997, Valtakunnan metsien… 2009). Valta pituus on osoittautunut varsin luotettavaksi lento- laserkeilauksella saatavaksi tunnukseksi (Næsset 2002, Järnstedt 2010, Nord-Larsen ja Riis-Nielsen 2010). Siksi vaihtoehtoisia malleja esitetään myös valtapituuden funktiona. Pahimmassa tapaukses- sa metsikön rakenteesta tunnetaan vain leimikon hakkuusopimuksen mukaiset tunnukset eli tukki- ja kuitupuun tilavuudet puulajeittain. Myös näiden tilavuustunnusten pohjalta laadittiin vaihtoehtoinen pituuden ennuste. Lopulliset mallit estimoitiin seka mallina R-ohjelmiston nlme-kirjaston lineaarisen sekamallin lme-funktiolla. 3 Tulokset 3.1 Ennustamismenetelmien vertailua Pituusmallien estimoimiseksi on käytettävissä eri- laisia vaihtoehtoisia lähestymistapoja. Vaihtoehtoi- silla menetelmillä A, B ja C estimoitujen mallien luotettavuuksia verrattiin aluksi keskenään. Koska mallien satunnaisosat ja virhetermit menetelmien välillä kuvasivat aivan eri asioita, ne jätettiin taulu- kosta pois. Vaikka menetelmät erosivat toisistaan, niin silti kiinteän osan estimoidut vakiot ja kertoimet poikkeavat yllättävän paljon toisistaan (taulukko 2). Esimerkiksi menetelmät A ja C ennustavat samoja parametreja logaritmisessa muodossa ja erona on mallin virheen minimointi puustotunnusten (mene- telmä A) tai puun pituuden (menetelmä C) mukaan. Menetelmän 1 logaritmimuunnoksilla aikaansaatu parametrien ja puustotunnusten välisten riippuvuuk- sien linearisointi paransi parametrien ennustemallia suhteessa Siipilehdon (1999) malliin (kuva 2). Silti samansuuntaiset trendit pituuden harhassa (havait- tu – ennuste) oli havaittavissa. Pituus yliarvioitiin yleisesti keskipituuden ollessa 13–20 metriä (kuva 2a). Keskipituuden ylittäessä 23 m puutason mallien harhat poikkesivat selvästi toisistaan siten, että me- netelmä 2 selvästi yliarvioi ja menetelmä 3 lievästi aliarvioi pituutta (kuva 2a). Puun suhteellisen koon mukaan harhat olivat toistensa kaltaisia, mutta niissä oli pienet tasoerot (kuva 2b). Molemmat puutason mallit (menetelmät B ja C) olivat parempia kuin pituuskäyrän parametrien ennustemallit (menetelmä A ja Siipilehto 1999). Molemmat puutason sekamallit antoivat melko hyvän tuloksen, mutta tässä tapauksessa suora pi- tuuden ennuste oli muunnosta parempi vaihtoehto (kuva 2). Linearisoidun mallin etuina oli kuitenkin metsikön sisäinen jäännöshajonnan homoskedasti- 222 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli suus ja ennen kaikkea satunnaistermien estimoin- ti voidaan tehdä lineaarisen ennustamisen teorian mukaisesti. Sen sijaan epälineaarisessa mallissa sa- tunnaisparametrien ennustaminen vaatii iteratiivisen ratkaisun (ks. Sirkiä ym. 2015). Siksi lineaarimuun- noksen (yhtälö 4) ennustaminen lineaarisella seka- mallilla valittiin lopullisten mallien lähtökohdaksi. Lisämuuttujilla ja muunnoksilla pyritään korjaa- maan kuvassa 2 havaittua harhaa ja aikaansaamaan aineiston vaihtelualueella mahdollisimman harhat- tomat ennusteet. Taulukko 2. Männiköiden pituuskäyrän ennusteet parametreille b0 ja b1 varttuneiden puustojen keskitunnuksilla DGM ja HGM. Mallin kiinteä osa menetelmällä A) sovitettujen pituuskäyrän parametrien ennustaminen, B) lineaari- sen muunnoksen ennustaminen (lme) ja C) pituuden ennustaminen epälineaarisesti (nlme, HDnaslund4 funktio, ks. Mehtätalo 2015). Siipilehto (1999) laati mallin suoraan b1 parametrille ja parametri b0 ratkaistiin siten, että käyrä kuli pisteen (DGM, HGM) kautta. Menetelmä: A B C Siipilehto (1999) b0 b1 b0 b1 b0 b1 b1 Vakio 1,559 0,533 1,422 0,449 0,995 0,658 Vakio 0,291 ln(DGM) 0,815 0,130 1,051 0,022 0,656 0,166 DGM 0,00134 ln(HGM) –1,048 –0,478 –1,277 –0,108 –0,679 –0,596 HGM –0,00634 –1,2 –1 –0,8 –0,6 –0,4 –0,2 0 0,2 0,4 0,6 0,8 1 1,2 5 10 15 20 25 Pi tu ud en h ar ha , m HGM,m a Menetelmä A Menetelmä B Menetelmä C Siipilehto 1999 –1 –0,8 –0,6 –0,4 –0,2 0 0,2 0,4 0,6 0,8 1 0,5 0,7 0,9 1,1 1,3 1,5 1,7 Pi tu ud en h ar ha , m d/D b Menetelmä A Menetelmä B Menetelmä C Siipilehto 1999 Kuva 2. Eri menetelmillä laadittujen mallien pituuden harhat keskipituu- den (HGM, kuva 2a) ja puun aseman (d/D, kuva 2b) suhteen. Menetelmän A parametrien logaritmimuunnoksilla ln(b0) ja ln(b1) saatiin selvä parannus Siipilehdon (1999) malliin, jossa b1 oli suoraan DGM ja HGM tunnusten line- aarinen funktio ja b0 ratkaistiin siten, että käyrä kuli pisteen (DGM, HGM) kautta. Molemmat puutason mallit B ja C olivat parempia kuin pituuskäy- rän parametrien ennustemallit. Myös C mallintaa logaritmisia parametreja (Mehtätalo 2015, HDnaslund4). 223 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… 3.2 Ennustemallit On hyvin ymmärrettävää, että keskiläpimitta ja keskipituus yhdessä selittävät hyvin pituuskäyrän parametreja, koska keskitunnukset kuvaavat läpi- mitta–pituus-koordinaatistossa käyrän jyrkkyyttä. Taimikoiden ja nuorten puustojen tunnuksilla ennus- tettaessa mallin tarkkuus kärsii melko selvästi, kun keskiläpimitta (D) ei ole käytettävissä selittävänä muuttujana (taulukko 3). Sen sijaan valtaläpimittaa (DDOM) harvoin tunnetaan ja se osoittautui huo- noksi selittäväksi muuttujaksi. Harhattomampi malli saatiin valtapituuden ja pohjapinta-alan funktiona (taulukko 5). Laadituissa Näslundin pituuskäyrän sekamalleissa (taulukot 3–7) metsikkötason vakion keskihajonta kasvoi selvästi, kun siirryttiin aritmeet- tisista keskitunnuksista painotettuihin keskitunnuk- siin ja siitä edelleen valtapituuteen. Tilavuustunnuk- silla ennustettaessa metsikkötason vakion hajonta oli odotetusti suurin. Mitä suurempi metsikkötason satunnaisvaikutus on, sitä merkittävämpää on mal- lin lokalisointi mittaushavaintojen avulla. Liitteessä esitetään esimerkki leimikkotunnuksilla ennustetun pituuskäyrän lokalisoimisesta (ks. taulukko 7). Jäännösvaihtelun keskihajonta kasvoi lämpö- summan pienetessä, mutta kun varianssifunktion selittävä muuttuja muotoiltiin 1000/ddy, tuli es- timoidusta potenssista positiivinen (0,77–1,09). Solakkuustermeille estimoitu potenssi oli noin –5 männylle, –5,8 kuuselle ja –4,5 koivulle merkiten voimakasta varianssin pienenemistä solakkuuden kasvaessa (taulukot 4 ja 5). Varianssifunktion esti- moidut potenssit muuttuivat melko vähän eri selit- täjien välillä, mutta melko paljon puulajien välillä (taulukot 3–7). Puuston kehitysvaihe vaikutti met- sikkökohtaisen jäännöshajontaan, mutta se kuvautui osittain solakkuuden kautta. 3.3 Mallin havainnollistaminen esimerkki­ metsikön avulla Männyn pituuskäyrän visuaaliseen tarkasteluun valittiin sama INKA–metsikkö (koe 1010, mittaus- kerta 1, lämpösumma 1267°Cvrk) kuin Siipilehdon (2011b, kuva 3) tutkimuksessa. Tähän 24 vuotiaa- seen VT männikköön ennustettiin pituuskäyrä vaih- toehtoisilla puustotunnuksilla, joita oli joko (N, D, H), (G, DGM, HGM) tai (G, HDOM). Metsikön koepuista oli maastossa mitattu läpimitta ja pituus. 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 16 18 20 Pi tu us , m Läpimitta, cm INKA 1010 D, H DGM, HGM DDOM, HDOM h(N, D, H) h(G, DGM, HGM) h(G, HDOM) Kuva 3. Pituuskäyrän ennusteita varttuneen männyn taimikkoon vaihtoehtoisilla puusto- tunnuksilla: h = f(N, D, H) (─ · ─); h = f(G, DGM, HGM) ( — ) ja h = f(G, HDOM) (- - -). Mitatut koepuut (u) ja lukupuut, joille on generoitu pituus metsikköön soviteun pituuskäyrän ja jäännösvaihtelun hajonnan avulla (◆). Lisäksi kuvassa on metsikköä kuvaavat puusto- tunnukset ja ±2 × sh kuvaamassa 95 % ennustetusta pituuden satunnaisvaihtelusta (····). 224 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli Taulukko 3. Yleiset Näslundin pituuskäyrän sovitukset männylle, kuuselle ja lehtipuustolle (ks. Kangas ja Maltamo 2002). Mallin satunnaistermit olivat metsikkötason i vakio (b0i), metsikkö- tason kerroin (b1i), metsikön mittauskertatason vakio (mitk) ja puutason jäännösvirhe (eij) ja taulukossa on niiden keskihajonnat s(). Termi s(eij) on varianssifunktion skaalausparametri ja varianssifunktion mukainen keskihajonta on siten syi = s(eij) (1000/DDY)p = 0,244 (1000/DDY)1,014. Parametri Mänty Kuusi Koivu b0 b0 b1 b0 b1 b0 b1 Estimaatti 1,181 0,247 1,706 0,327 1,014 0,238 Keskivirhe 0,015 0,002 0,026 0,003 0,044 0,008 s(b0i) 0,299 0,339 0,296 s(b1i) 0,046 0,038 0,059 corr(b0, b1) –0,232 –0,683 –0,778 s(mitk) 0,184 0,057 0,067 s(eij) 0,244 0,262 0,203 varianssifunktio (1000/DDY) 1,014 0,779 1,026 Taulukko 4. Pituuskäyrän ennusteet parametreille b0 ja b1 nuorten puustojen tunnuksilla keskiläpimitta (D) tuntien tai ilman sitä. Kaikki selittävät muuttujat olivat tilastollisesti erittäin merkitseviä. Varianssifunktio laadittiin suhteessa solakkuuteen ln(H/D+2), kun sekä D että H tunnettiin ja muulloin lämpösummaan (1000/DDY). Näiden estimoidut potenssit annetaan kah- della viimeisellä rivillä. Parametri Mänty Kuusi Koivu Estim. Estim. Estim. Estim. Estim. Estim. b0 Vakio 1,481 5,120 –0,585 6,519 0,824 1,886 ln(DDY) –0,186 –0,345 –1,045 √T 0,067 1/√T –2,111 1,397 6,136 ln(T) 0,557 ln(N) –0,126 –0,187 –0,136 –0,110 ln(D) 2,279 1,741 0,837 H 0,031 0,067 0,078 H2 0,00145 H/D 1,744 ln(H–1) –2,297 –0,627 –1,117 –0,477 ln(H–2) –1,036 b1 Vakio 0,233 0,348 0,491 0,451 0,470 0,485 1/√T 0,287 0,366 0,369 ln(N) 0,010 ln(D) 0,041 ln(H–1) –0,091 –0,069 –0,106 –0,111 ln(H–2) –0,070 –0,070 s(b0i) 0,153 0,191 0,244 0,307 0,084 0,220 s(b1i) 0,014 0,016 0,019 0,022 0,013 0,022 corr(b0, b1) –0,908 –0,786 –0,964 –0,898 –0,959 –0,833 s(eij) 0,307 0,248 0,309 0,260 0,343 0,205 varianssifunktio ln(H/D+2) –5,001 –5,831 –4,747 (1000/DDY) 0,992 0,800 0,944 225 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… Taulukko 5. Pituusmallit ja varianssifunktio varttuneiden puustojen puustotunnuksilla ennustettuna. Parametri Mänty Kuusi Koivu Estim. Keskivirhe Estim. Keskivirhe Estim. Keskivirhe b0 Vakio 3,394 0,211 3,320 0,536 0,853 0,074 ln(DDY) –0,296 0,027 –0,365 0,068 ln(G+1) –0,253 0,018 ln(DGM) 0,769 0,061 0,849 0,173 0,525 0,093 HGM 0,029 0,004 0,057 0,007 0,058 0,008 ln(HGM–1) –0,838 0,060 –0,946 0,156 –0,819 0,106 b1 Vakio 0,431 0,006 0,584 0,012 0,504 0,012 ln(G+1) 0,017 0,001 ln(DGM) 0,034 0,005 0,047 0,011 0,046 0,010 ln(HGM–1) –0,136 0,005 –0,147 0,010 –0,161 0,011 s(b0i) 0,182 0,252 0,080 s(b1i) 0,013 0,015 0,0057 corr(b0, b1) –0,914 –0,956 –0,872 s(eij) 0,275 0,291 0,206 Varianssifunkio ln(HGM/DGM+2) –5,349 –5,842 –4,476 Taulukko 6. Pituusmallit valtapituuden avulla ennustettuna. Parametri Mänty Kuusi Koivu Estim. Keskivirhe Estim. Keskivirhe Estim. Keskivirhe b0 Vakio 7,136 0,321 9,833 0,759 0,206 0,166 ln(DDY) –0,686 0,045 –1,185 0,106 –0,708 0,128 ln(G+1) –0,273 0,019 HDOM2 0,00139 8,2 10–5 0,00176 1,3 10–4 9,15 10–4 2,0 10–4 ln(HDOM–1) –0,329 0,029 –0,188 0,103 ln(HDOM–4) 0,098 0,044 b1 Vakio 0,530 0,0053 0,710 0,013 0,474 0,014 ln(G+1) 0,0136 0,0015 ln(HDOM–1) –0,128 0,0024 –0,133 0,0046 ln(HDOM–4) –0,104 0,006 s(b0i) 0,235 0,275 0,177 s(b1i) 0,015 0,016 0,018 corr(b0, b1) –0,765 –0,797 –0,639 s(eij) 0,249 0,259 0,204 Varianssifunkio (1000/DDY) 0,994 0,813 1,056 Lukupuista oli mitattu vain läpimitta ja pituus ge- neroitiin metsikköön sovitetun Näslundin käyrän ja metsikkö/mittauskerta kohtaisen, aineistosta es- timoidun hajontatermin avulla. Koepuut ja lukupuut on kuvassa erotettu toisistaan. Koko puujoukosta lasketut metsikön puustotunnukset olivat: D 9,2 cm, H 8,4 m, DGM 10,4 cm, HGM 9,0 m, DDOM 14,3 cm, HDOM 9,9 m ja tiheystunnukset N 2 577 ha–1 ja G 18,4 m2ha–1. Metsikön keskitunnuksilla (D, H) ja (DGM, HGM) 226 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli ennustetut pituuskäyrät sopivat erinomaisesti ha- vaintoaineistoon (kuva 3). Molemmat käyrät näyt- tivät kulkevan lähes selittäjänä olleen pisteen kautta. Valtapituudella ennustettu pituuskäyrä puolestaan aliarvioi pituutta tässä metsikössä ja kulki pisteen (DDOM, HDOM) alapuolelta. (Kuvan ennusteissa on mukana harhankorjaus). Metsikön puiden pituu- det pysyivät ennustetun satunnaisvaihtelun rajoissa, joka on kuvattu pisteviivoilla (2 × sh eli 95 %:n rajat). Tässä kuvassa vaihteluväli edusti HGM/DGM suh- teen avulla ennustettua jäännöshajontaa (sy(HGM, DGM) = 0,240). Aritmeettisilla tunnuksilla tai läm- pösumman avulla ennustetut vaihteluvälit olisivat olleet melko samanlaiset, koska sy(H, D) = 0,228 ja sy(DDY) = 0,243. 3.4 Näslundin pituuskäyrän lokalisointi Koska ennustemallit laadittiin lineaarisina seka- malleina, ne olivat lokalisoitavissa (kalibroitavissa) yhden tai useamman koepuumittauksen avulla line- aarisen ennustamisen teorian mukaan (esim. Lappi 1991). Tämän tutkimuksen liitteessä on esimerkki, jossa yksi malli lokalisoidaan viiden ensimmäi- sen koepuumittaushavainnon avulla. Lokalisoin- nissa käytetään ns. parasta lineaarista harhatonta ennustinta (BLUP) mallin satunnaisparametreille. Liitteen esimerkissä käytettiin epävarminta mallia, jossa Näslundin pituuskäyrä ennustettiin puuston tilavuustunnusten avulla ilman yhtään keskitun- nusta (metsänhakkuusopimuksen mukaiset lähtö- tiedot, taulukko 7). On selvää, ettei puuston tila- vuus anna yksiselitteistä kuvaa puuston rakenteesta. Kuitupuun ja tukkipuun osuudet antoivat osviittaa puuston kehitysvaiheesta, mutta pituuskäyrän en- nustetarkkuudessa oli suurta vaihtelua metsiköiden välillä. Esimerkkitapauksessa kiinteän osan ennuste yliarvioi puiden pituuksia selvästi. Esitetyillä loka- lisoinnin vaihtoehdoilla, eli 1, 3 tai 5 ensimmäistä mittaushavaintoa, saatiin pituuskäyrää oleellisesti korjattua, vaikka ensimmäinen havainto ei vielä riittänyt korjaamaan pituuskäyrän muotoa (ks. lii- te1 ja liitekuva 1). 3.5 Mallien testausta Tarkastellaan lähemmin yleisimmin sovellettavaa mallia, jossa tunnetut puustotunnukset ovat tyypil- lisiä varttuneen puuston tunnuksia, eli G, DGM ja HGM. Residuaalikuvat on laskettu mallinnusaineis- Taulukko 7. Pituusmallit leimikon tukkilavuuden ja kuitupuutilavuuden avulla ennustettuna. Kaupallinen tilavuus (Vkaup) on tukki- ja kuitupuutilavuuden summa. Parametri Mänty Kuusi Koivu Estim. Keskivirhe Estim. Keskivirhe Estim. Keskivirhe b0 Vakio 3,128 0,065 3,011 0,144 1,400 0,161 (DDY/1000) –0,537 0,058 –0,963 0,103 –0,666 0,158 ln(Tukki+2) –0,041 0,005 0,161 0,013 ln(Kuitu+2) –0,414 0,010 –0,212 0,026 –0,346 0,101 ln(Vkaup+2) 0,409 0,097 b1 Vakio 0,232 0,003 0,387 0, 009 0,343 0,010 ln(Tukki+2) –0,025 0,0005 –0,021 0,001 ln(Kuitu+2) 0,023 0,001 0,006 0,002 0,023 0,006 ln(Vkaup+2) –0,054 0,006 s(b0i) 0,245 0,295 0,209 s(b1i) 0,022 0,022 0,025 Corr(b0, b1) –0,621 –0,790 –0,673 s(eij) 0,257 0,261 0,209 Varianssifunktio (1000/DDY) 0,914 0,800 0,790 227 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… tosta. Lämpösumman suhteen ei ollut havaittavissa minkäälaista trendiä millään puulajilla (kuva 4 A). Suurimmillaan tavattiin pituuden 31 cm:n aliarvio koivulla lämpösumman ollessa noin 1100 °Cvrk. Pituusmalleissa ei ollut myöskään harhan trendiä keskipituuden suhteen (kuva 4B). Kuvassa 2 ha- vaittua harhaa saatiin korjattua pienillä muutoksil- la, kuten lisäämällä HGM sellaisenaan selittämään parametria b0 ja toisaalta muunnos ln(HGM – 1) oikaisi residuaalia yllättävän paljon verrattuna aluksi kokeiltuun muuttujaan ln(HGM). Silti aivan kook- kaimmissa kuusikoissa pituuden 38 cm yliarvio oli huomionarvoista. Suhteellinen virhe oli silti vain noin 1,4 %. Havaintoja kyseisessä luokassa oli vain 78, kun niitä tyypillisesti oli yli 300. Puun aseman (d/D) suhteen on havaittavissa pituuden yliarviota metsikön pienimmissä puissa ja aliarviota metsikön keskikokoisissa puissa (kuva 4C). Suurimmillaan aliarvio oli männyllä 13 cm, kuusella 8 cm ja koi- vulla 11 cm. Männyllä kookkaimpien puiden (d/D = 1,7) pituus hieman yliarvioitiin, mutta kuusella ja koivulla nämä ennusteet kääntyivät edellisen luo- –0,6 –0,4 –0,2 0,0 0,2 0,4 0,6 850 900 950 1000 1050 1100 1150 1200 1250 1300 Pi tu ud en h ar ha , m Lämpösumma A Mänty Kuusi Koivu –0,6 –0,4 –0,2 0,0 0,2 0,4 0,6 5 7 9 11 13 15 17 19 21 23 25 27 Pi tu ud en h ar ha , m HGM, m B –0,6 –0,4 –0,2 0,0 0,2 0,4 0,6 0.5 0,7 0,9 1,1 1,3 1,5 1,7 Pi tu ud en h ar ha , m d/D C Kuva 4. Puun pituuden harhat puulajeit- tain A: lämpösumman, B: mediaanipituu- den (HGM) ja C: puun suhteellisen koon (d/D) suhteen. Pituuskäyrät ennustettiin varttuneiden puustojen tunnuksilla. Y- akseli on skaalattu puoleen keskimää- räisestä 1,2 m hajonnasta. 228 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli kan 27–46 cm yliarvioista lähes harhattomiksi (kuva 4C). Kaiken kaikkiaan mallien harhat olivat pieniä verrattuna virheen keskihajontaan eri luokissa. Esimerkiksi lämpösummaluokissa keskihajonta oli 97–150 cm, d/D-luokissa keskihajonta oli 88–185 cm ja HGM-luokissa 52–220 cm. Keskimääräinen virheen keskihajonta oli noin 1,2 m ja harhakuvien y-akseli skaalattiin puoleen siitä. Mainittakoon, että keskimääräiset harhankorjaukset koko aineistossa olivat männyllä ja kuusella 12 cm ja koivulla 17 cm ja suurimmillaan ne olivat 39, 18 ja 29 cm. 4 Tulosten tarkastelua Näslundin (1936) pituuskäyrä on melko yksinkertai- nen, koska siinä on vain kaksi estimoitavaa paramet- ria. Kaksiparametrisista pituuskäyristä se on ollut yksi parhaimmista (ks. Fang ja Bailey 1998, Leduc ja Goelz 2009). Toisinaan se on sopinut eri puula- jien aineistoihin jopa paremmin kuin suositut kol- men parametrin funktiot (ks. Mehtätalo ym. 2015). Kun malli estimoitiin pituuden tai sen muunnoksen virheen minimoimiseksi, niin pituuden ennuste oli luotettavampi kuin tapauksissa, joissa minimoitiin etukäteen sovitettujen parametrien ennustevirhettä (Siipilehto 1999, 2011a, 2011b) ilman puiden pi- tuushavaintoja (ks. kuva 2). Tulos on analoginen jakaumaparemetrien ennusteiden kanssa. Tyypilli- sesti jakaumaparametrien ennustemallit laaditaan etukäteen sovitetuille parametreille, mutta Cao (2004) kehitti menetelmän, jossa malli sovitettiin maksimum likelihood -funktion osana suoraan puu- tason havaintoihin ja mallin luotettavuus samalla parani (ks. Siipilehto 2009). Jos mallitetaan pelkkiä pituuskäyrän parametreja ei estimoiduista virheter- meistä voi tehdä johtopäätöksiä itse puun pituuden tarkkuudesta. Näslundin pituuskäyrää sovelletaan sekä toisen että kolmannen asteen yhtälönä. Kuten Petterson (1955), Vestjordet (1972) ja Siipilehto (1999), myös tässä tutkimuksessa käytettiin kolmannen asteen yh- tälöä kuuselle, kun taas männyllä ja koivulla käy- tettiin Näslundin (1936) alkuperäistä toisen asteen yhtälöä. Kuusi poikkeaa varjoa sietävänä puulajina männystä ja koivusta. Kuusikoille on ominaista laaja läpimitta- ja pituusvaihtelu metsikön sisällä, joten kolmannen asteen yhtälö loi toisen asteen yhtälöä paremmat edellytykset sigmoidin pituuskäyrän so- vittamiseksi ja ennustamiseksi. Kolmannen asteen yhtälöä on käytetty myös Tanskassa varjoa sietä- vän pyökin pituuskäyrän tasoittamiseksi (Johannsen 2002, Nord-Larsen ja Cao 2006). Fitje ja Vestjordet (1977), Staudhammer ja LeMay (2000) ja Leduc ja Goelz (2009) kokeilivat Näslundin kolmannen asteen yhtälöä yhtenä vaihtoehtona, mutta päätyivät lopulta soveltamaan muita pituusmalleja. Schmidt ym. (2011) käyttivät kolmannen asteen yhtälöä männylle, kun taas Kinnunen ym. (2007) estimoi- vat toisen asteen yhtälön kuuselle ja muut puula- jit olivat yhteisen pituusmallin dummy-muuttujia. KPL-laskentaohjelmassa Näslundin pituuskäyrän potenssi on parametrisoitu eli sitä voidaan muut- taa, mutta potenssin valintaa ei ole toistaiseksi oh- jeistettu ja sen oletusarvo on puulajista riippumatta 2 (Heinonen 1994, s. 31). KPL-laskentaohjelman ohjeita ollaan päivittämässä ja samalla Näslundin yhtälön parametrisointi saa uudet ohjeet (Jaakko Heinonen suullisesti). Tässä tutkimuksessa laaditut pituusmallit olivat varsin luotettavia johtuen mm. siitä, että puustotun- nukset olivat tarkkaan mitatuista puista laskettuja eli aineiston puustotunnuksissa ei ollut tyypillistä ar- viointi virhettä. Ennusteen luotettavuuden kannalta onkin ensiarvoisen tärkeää, että malliin syötettävät puustotunnukset, kuten keskiläpimita ja keskipi- tuus ovat mahdollisimman tarkkoja. Silmävaraista arviointia tehdään laajemmassa mittakaavassa enää VMI:n yhteydessä, mutta sen tarkkuus vaikuttaa yhä kilpailukykyiseltä kuvatulkintamenetelmiin verrattuna (ks. Uuttera ym. 2002, Haara ja Korho- nen 2004, Uuttera ym. 2006). Laseraineistoilla on päästy parhaimmillaan alle 10 % keskivirheeseen metsikkötason keski- ja valtapituudessa tai keskilä- pimitassa (Næsset 2002, Suvanto ym. 2005, Uuttera ym. 2006, Packalén ja Maltamo 2007), mutta koea- latasolla tulokset eivät ole yhtä hyviä. Esimerkiksi Holopainen ym. (2009) ja Järnstedt (2010) vertasivat korkean resoluution ilmakuvia (0,5 m resoluutio) laserkeilaukseen (1–3 m resoluutio) ja saivat keski- läpimitan ja keskipituuden keskivirheiksi 25–34 % ilmakuvilta ja 20–29 % laseraineistolla koealatasol- la tarkasteltuna. Edellä mainituista Næsset (2002), Suvanto ym. (2005) ja Uuttera ym. (2006) käyttivät puustotunnusten regressiomallitusta laseraineiston 229 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… korkeusjakaumatunnuksista ja muissa sovellettiin laseraineiston aluepohjaista tulkintaa k-nn mene- telmällä. Kaiken kaikkiaan metsätalouden suunnit- telujärjestelmissä mallivirheiden (puuston luominen jakaumamallien ja läpimitta-pituus riippuvuuden avulla sekä kasvun ennustaminen) osuus on hyvin pieni verrattuna puustotunnusten arviointivirheeseen (ks. Ojansuu ym. 2002, Holopainen ym. 2010). Laaditut pituuskäyrän ennusteet kulkivat lähes mallien syöttötietona annetun pisteen kautta, oli pis- te sitten (D, H) tai (DGM, HGM). Siten lähtötiedon (mallien selittäjät) ja lopputuloksen (malleilla luodut kuvauspuut) välinen yhteensopivuus vaikuttaa hy- vältä. Mallit laadittiin lineaarisina sekamalleina, jotta ne voidaan helposti lokalisoida, eli ennustaa sekamallin satunnaisparametrit koepuumittausten avulla. Alustavan tarkastelun perusteella nyt laa- dittujen mallien kiinteä osa antaa hyvän tarkkuu- den sellaisenaan, joten lokalisoinnille ei ole yleistä tarvetta. Poikkeuksena leimikon tilavuustunnuk- silla ennustaminen (taulukko 7), jonka yhteydessä esitettiin esimerkki mallin lokalisoinnista. Samoin lokalisontia tarvitaan, kun tehdään yleinen malli pelkästään puun läpimitan funktiona (taulukko 3), kuten Kangas ja Maltamo (2002). Sekamallit loka- lisoituvat varsin tehokkaasti ja yksikin mittausha- vainto voi selvästi parantaa ennustetta (esim. Lappi 1991, Eerikäinen ym. 2002, Mehtätalo 2005). Jos havainto poikkeaa huomattavasti mallin odotusar- vosta, niin yksi havainto ei riitä korjaamaan ennus- tetta täysimääräisesti. Otetaan esimerkiksi INKA metsikön 1215 toinen mittauskerta, jossa mediaani DGM = 16,0 cm, HGM = 16,7 m (G = 12,5 m2ha–1, DDY 1086 °Cvrk) ja yleisen mallin (taulukko 3) pi- tuuden odotusarvo 11 m, kun läpimitta on 16 cm. Kun malli kalibroitiin mediaanipuun dimensioilla (ks. Kangas ja Maltamo 2002), saatiin mediaaniläpi- mittaa vastaavaksi pituudeksi 16,0 m. Kun metsikön puustotunnukset syötettiin regressiomalliin (tauluk- ko 5), niin mediaanipuun pituudeksi saatiin 16,7 m, eli tässä tapauksessa sama kuin mitattu HGM. Toisaalta, jos yksittäinen mittaushavainto on täysin sattumanvarainen, voi ennusteen lokalisoituminen olla melko vaatimatonta, kuten tämän tutkimuksen liitteessä esitetyssä esimerkkitapauksessa. Kun mit- taushavaintoja lisättiin kolmeen ja lopulta viiteen, niin lokalisoitu pituuskäyrä sopi aina vaan parem- min metsikön kaikkiin läpimitta- ja pituusmittaus- havaintoihin. Vaihtoehtoisia läpimitta-pituusmalleja Suomessa ovat esittäneet Lappi (1991, 1997), Hökkä (1997), Mehtätalo (2004, 2005) ja Eerikäinen (2009) ja useimmiten näissä on esitetty lokalisointi lineaari- sen ennustamisen sovelluksena mitattujen koepui- den avulla. Alustava tarkastelu osoitti, että epäline- aarinen funktio puun pituudelle, jossa mallitettiin logaritmisia parametreja (parametri on aina positii- vinen) selittäjien lineaarisella funktiolla (ks. Mehtä- talo 2015), tarjosi varteenotettavan vaihtoehdon. Jos oltaisiin kiinnostuneita pelkästään mallin kiinteästä osasta, se voisi olla paras vaihtoehto. Epälineaarisen mallin lokalisointi on kuitenkin hankalampaa, mutta se voidaan tehdä iteratiivisesti Taylorin sarjaan pe- rustuen (ks. Sirkiä ym. 2014). Tässä tutkimuksessa laadituilla vaihtoehtoisilla malleilla tulisi korvata Siipilehdon (1999) Näslundin pituuskäyrän parametreja ennustavat mallit mahdol- lisissa suunnittelujärjestelmissä tai tutkimusaineis- toja malleilla täydennettäessä. Uusilla malleilla voi- daan välttää Siipilehdon (2011b) havaitsema pituu- den harha, joka johtui sekä Näslundin parametrien ennustamisesta ilman logaritmimuunnoksia Siipilehdon (1999) mallilla, että kyseisen mallin soveltamisesta laadinta-aineiston ulkopuolella. Tämä harha oli havaittavissa erityisesti nuorissa metsiköissä sekä varttuneiden metsiköiden vallitun latvuskerroksen aluspuissa (Siipilehto 2011b, kuva 4, s 30). Kuusisto ja Kangas (2008) havaitsi Siipilehdon (1999) kuusen pituusmallin aiheuttavan tilavuustunnusten kasvavaa aliarviota keskipituuden pienetessä. Nyt laadittu kuusen pituusmalli oli lähes harhaton keskipituuden suhteen (kuva 4b) ja pie- nimmissä kuusikoissa (HGM 5–7 m) aliarvio oli vain1–2 %. Kun pienimmät puustot ennustettiin nuorten puustojen tunnuksilla (N, D, H), niiden tark- kuus edelleen parani (0,3–0,7 %). Taulukoiden 4–7 malliperheestä voidaan valita kulloisenkin käyttäjän lähtötietoihin parhaiten sopiva vaihtoehto. Pituuskäyrän lisäksi malleilla voitiin ennustaa pituuden metsikkökohtaista keskihajontaa vari- anssifunktioilla. Aikanaan Siipilehto (2000) esitti yksinkertaiset mallit keskihajonnan ennustamiseksi perustuen solakkuuteen (HGM/DGM). Myös tässä tutkimuksessa puuston keskimääräinen solakkuus osoittautui tärkeimmäksi jäännöshajontaa selittäväk- si tunnukseksi. Jos keskitunnuksia ei voitu olettaa 230 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli tunnetuksi, varianssifunktio kuvattiin lämpösumman avulla. Jäännösvirheen keskihajonta kasvoi selvästi etelästä pohjoiseen. Esimerkiksi männyllä lineaari- muunnoksen y jäännösvirheen keskihajonnaksi saa- tiin varianssifunktiolla 0,22, kun lämpösumma oli 1200 °Cvrk ja 0,27, kun lämpösumma oli 800 °Cvrk. Toisaalta, jos katsotaan puuston solakkuutta medi- aaniläpimitan DGM ollessa 20 cm ja pituuden HGM vaihdellessa 12 m:stä 24 m:iin, niin näitä vastaavat keskihajonnat männyllä olivat varianssifunktion mukaan 0,31 ja 0,18. Harhankorjauksen lisäksi hajontatieto on oleel- linen silloin, kun halutaan generoida satunnaista vaihtelua pituuden odotusarvon ympärille eli läpi mitan suhteen ehdollista pituusvaihtelua. Pituuden satunnaistamisella voidaan saavuttaa joitakin sel- keitä hyötyjä. Satunnaistamisen avulla pituuden reunajakauma vastaa hyvin luonnossa havaittavaa pituuden reunajakaumaa, kun taas pituuden odotus- arvo supistaa sitä voimakkaasti (Siipilehto 2000). Siipilehto (2001) osoitti, että näin käy etenkin van- hoissa männiköissä, joissa läpimittaluokan sisäinen pituusvaihtelu on selvästi suurempaa kuin läpimit- taluokkien välinen pituusvaihtelu (ks. kuva 5). Lä- pimittaluokan sisäisestä pituusvaihtelusta voi olla hyötyä myös harvennusten generoimisessa, koska esimerkiksi alaharvennus perustuu paremminkin puiden pituuksien välisiin suhteisiin kuin puun lä- pimittaan (Hafley ja Buford 1985). Myöskin realisti- sen testiaineiston luomiseksi satunnaistettua pituutta on pidetty parempana vaihtoehtona kuin pituuden odotusarvo (Siipilehto 2011a, 2011b). Metsikön kokonaistilavuuteen ja puutavaralajien tilavuus- tunnuksiin sillä ei sen sijaan ole juuri merkitystä (Siipilehto 2000). Näslundin pituuskäyrä on suosittu koepuiden pituuden tasoittamiseksi ja yleistämiseksi puil- le, joiden pituutta ei tunneta. Toisaalta Näslundin pituusmallia on käytetty apumallina esimerkiksi tilavuuden arvioimiseksi, eikä pituusmallin hyvyyttä ole näissä tutkimuksissa erikseen tarkasteltu (esim. Siipilehto 1999, Kangas ja Maltamo 2002, Kinnunen ym. 2007, Vastaranta ym. 2010, Siipilehto 2011a). Siten Näslundin pituuskäyrän ennustamisesta ja en- nusteen tarkkuudesta on toistaiseksi melko vähän käytännön esimerkkejä. Nyt laaditut mallit eivät olleet täysin harhattomia, mutta tyypillisillä puus- totunnuksilla saatiin hyviä tuloksia. Tämä tarkoittaa mm. sitä, että N, D ja H olivat hyviä selittäjiä nuo- rille puustoille, kun taas G, DGM ja HGM toimivat hyvin varttuneille puustoille. Myös valtapituuden ja lämpösumman (männyllä myös pohjapinta-alan) funktiona pituuskäyrä saatiin kohtalaisen luotetta- vasti. Viimeksi mainittuja malleja voitaneen tule- vaisuudessa soveltaa laseraineiston yhteydessä, kos- ka valtapituus saadaan luotettavasti laserpisteiden korkeusjakaumasta (Næsset 2002, Nord-Larsen ja Riis-Nielsen 2010). Toistaiseksi suomalaiset sovel- lukset valtapituuden ennustamiseksi laseraineistosta puuttuvat. Heikoimmillaan metsikön kuvaus voi si- sältää vain puulajeittaiset puutavaralajien tilavuudet eli metsänhakkuusopimuksen mukaiset tunnukset. Tällaisessa tilanteessa on aika mahdotonta saada luotettavaa puuston rakenteen kuvausta. Siksi tähän tilanteeseen sopiva sovellus olisi hakkuukoneen kor- juuaikainen mallien kalibrointi. Tällaista tapausta si- muloitiin liitteessä ja mallin lokalisointi esim. viiden ensimmäisen puun jälkeen tuotti varsin luotettavan pituuden ennusteen. 0 5 10 15 20 25 30 0 10 20 30 40 P itu us , m Läpimitta, cm Kuva 5. Esimerkki läpimitta-pituus -riippuvuudesta vanhassa männikössä (T = 70,  DGM = 27 cm). Läpimitta havaintojen (x) vaihtelualueella männyn ennustettu pituus (—) vaihteli välillä 22,5–26 m. Ennustettu satunnaisvaihtelu (- - -) lisäsi vaihteluvälin 19–28,5 metriin eli lähelle havaittua vaihtelua. Kuvassa vaihteluväli on pituuden odotusarvo ±2 × sh eli se edustaa 95 % satunnaisvaihtelusta. 231 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… Kirjallisuus Arabatzis, A.A. & Burkhart, H.E. 1992. An evaluation of sampling methods and model forms for estimating height-diameter relationships in loblolly pine planta- tions. Forest Science 38: 192–198. Cao, Q.V. 2004. Predicting parameters of a Weibull func- tion for modeling diameter distribution. Forest Science 50(5): 682–685. Curtis, R.O. 1967. Height-diameter and height-diameter- age equations for second-growth Douglas-fir. Forest Science 13(4):265–357. Eerikäinen, K. 2009. A multivariate linear mixed-effects model for the generalization of sample tree heights and crown ratios in the Finnish National Forest Inventory. Forest Science 55(6): 480–493. — & Korhonen, K.T. 2001. Puutunnusmallit. Julkaisussa: Maltamo, M. & Laukkanen, J. (toim.). Metsää kuvaa- vat malli. Silva Carelica 36. Joensuun yliopisto 239 s. — , Mabvurira, D., Nshubemuki, L., Saramäki, J. 2002. A calibrateable site index model for Pinus kesiya planta- tion in southeastern Africa. Canadian Journal of Forest Research 32: 1916–1928. Elfving, B. & Kiviste, A. 1997. Constructing of site index equations for Pinus sylvestris L. Using permanent plot data in Sweden. Forest Ecology and Management 98: 125–134. Fahlvik, N., Elfving, B. & Wikström, P. 2014. Evaluation of growth models used in the Swedish forest planning system Heureka. Silva Fennica 48(2). 17 s. Article id 1013. Fang, Z. & Bailey, R.L. 1998. Height‐diameter models for tropical forests on Hainan Island in southern China. Forest Ecology and Management 110(1–3): 315–327. Fitje, A. & Vestjordet, E. 1977. Bestandshøydekurver og nya høydeklasser for gran. Meddelelser fra Norsk In- stitutt for Skogforskning 34(2). 68 s. Gómez-García, E., Diéguez-Aranda, U., Castedo-Dorado F. & Crecente-Campo, F. 2014. A comparison of model forms for the development of height-diameter rela- tionships in even-aged stands. Forest Science 60(3): 560–568. Gustavsen, H.G., Roiko-Jokela, P. & Varmola, M. 1988. Kivennäismaiden talousmetsien pysyvät (INKA ja TINKA) kokeet. Suunnitelmat, mittausmenetelmät ja aineistojen rakenteet. Metsäntutkimuslaitoksen tiedon- antoja 292. 212 s. Haara, A. & Korhonen. K.T. 2004. Kuvioittaisen arvioin- nin luotettavuus. Metsätieteen aikakauskirja 4/2004: 489–508. Hafley, W.L. & Buford, M.A. 1985. A bivariate model for growth and yield prediction. Forest Science 31(1): 237–247. Heinonen, J. 1994. Koealojen puu- ja puustotunnusten laskentaohjelma KPL. Käyttöohje. Metsäntutkimus- laitoksen tiedonantoja 504. 80 s. Holopainen, M., Tuominen, S., Karjalainen, M., Hyyppä, J., Vastaranta, M., Hujala, T. & Tokola, T. 2009. Kor- kearesoluutioisten E-SAR-tutkakuvien tarkkuus puus- totunnusten koealatason estimoinnissa. Metsätieteen aikakauskirja 4/2009: 309–323. — , Vastaranta, M., Rasinmäki, J., Kalliovirta, J., Mäki- nen, A., Haapanen, R., Melkas, T., Yu, X. & Hyyppä, J. 2010. Uncertainty in timber assortment estimates predicted from forest inventory data. European Journal of Forest Research 129: 1131–1142. Hökkä, H. 1997. Height-diameter curves with random intercepts and slopes for trees growing on drained peatlands. Forest Ecology and Management 97: 63–72. Johannsen, V.K. 2002. Selection of diameter-height curves for even-aged oak stands in Denmark, Dynamic growth models for Danish forest tree species, Working paper 16, Danish Forest and Landscape Research Institute. Järnstedt, J. 2010. Erittäin korkean resoluution ilmaku- vista tuotettu pintamalli puustotunnusten estimoinnis- sa. Metsänarvioimistieteen ja -teknologian pro gradu -tutkinto. Helsingin yliopisto. 52 s. Kangas, A. & Maltamo, M. 2002. Anticipating the varian- ce of predicted stand volume and timber assortments with respect to stand characteristics and field measu- rements. Silva Fennica 36(4): 799–811. — , Mäkinen, H. & Lyhykäinen, H.T. 2010. Value of qual- ity information of Scots pine stands in timber bidding. Canadian Journal of Forest Research 40: 1781–1790. — , Hurttala, H., Mäkinen, H. & Lappi, J. 2012. Esti- mating the value of wood quality information in con- strained optimization. Canadian Journal of Forest Research 42: 1347–1358. Kinnunen, J., Maltamo, M & Päivinen, R. 2007. Standing volume estimates of forests in Russia: how accurate is the published data? Forestry 80(1): 53–64. Kuusela, K. & Salminen, S. 1969. The 5th National Forest Invetory in Finland. General design, instruction for field work and data processing. Metsäntutkimuslai- toksen Julkaisuja 69(4). 232 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli Kuusisto, L. & Kangas, A. 2008. Harhakomponentit ku- vioittaisen arvioinnin puuston tilavuuden laskentaket- jussa. Metsätieteen aikakauskirja 3/2008: 177–190. Lappi, J. 1991. Calibration of height and volume equa- tions with random parameters. Forest Science 37(3): 781−801. — 1993. Metsäbiometrian menetelmiä. Silva Carelia 24. Joensuun yliopisto. 182 s. — 1997. A longitudinal analysis of height/diameter cur- ves. Forest Science 43(4): 555–570. Leduc, D. & Goelz, J. 2009. A height-diameter curve for lonfleaf pine plantations in the Gulf Coastal Plain. Southern journal of Applied Forestry 33(4): 164–170. Maltamo, M., Eerikäinen, K., Packalén, P. & Hyyppä, J. 2006. Estimation of stem volume using laser scanning- based canopy height metrics Forestry 79(2): 217–229. Mehtätalo, L. 2004. A longitudinal height-diameter model for Norway spruce in Finland. Canadian Journal of Forest Research 34: 131–140. — 2005. Height-Diameter (H-D) models for Scots pine and birch in Finland. Silva Fennica 39(1): 55–66. — 2015. Package lmfor. Saatavilla: https://cran.r-project. org/web/packages/lmfor/lmfor.pdf. — , de-Miguel, S. & Gregoire, T.G. 2015. Modeling height-diameter curves for prediction. Canadian Jour- nal of Forest Research 45(7): 826–837. Mielikäinen, K. 1980. Structure and development of mi- xed pine and birch stands. Communicationes Instituti Forestalis Fenniae 99(3). 82 s. — 1985. Effect of an admixture of birch on the structure and development of Norway spruce stands. Commu- nicationes Instituti Forestalis Fenniae 133. 73 s. Muinonen E., Anttila P., Heinonen J. & Mustonen J. 2013. Estimating the bioenergy potential of forest chips from final fellings in Central Finland based on biomass maps and spatially explicit constraints. Silva Fennica 47(4), article id 1022. 22 s. Mustonen, J., Packalén, P & Kangas, A. 2008. Automatic segmentation of forest stands using a canopy height model and aerial photography. Scandinavian Journal of Forest Research 23: 534–545. Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two- stage procedure and field data. Remote Sensing of Environment 80: 88–99. Newton, P.F. & Amponsah, I.G. 2007. Comparative eva- luation of five height-diameter models developed for black spruce and jack pine stand-types on terms of goodness-of-fit, lack-of-fit and predictive ability. Fo- rest Ecology and Management 247: 149–166. Nord-Larsen, T. & Cao, Q.V. 2006. A diameter distri- bution model for even-aged beech in Denmark. Forest Ecology and Management 231(1–3): 218–225. — & Riis-Nielsen, T. 2010.Developing an airborne laser scanning dominant height model from a countrywide scanning survey an national forest inventory data. Scandinavian Journal of Forest Research 25: 262–272. Näslund, M. 1936. Skogsförsöksanstaltens gallrings- försök i tallskog. Meddelanden från Statens Skogs- försöksanstalt 29. 169 s. Ojansuu, R., Halinen, M. & Härkönen, K. 2002. Metsä- talouden suunnittelujärjestelmän virhelähteet männyn ensiharvennuskypsyyden määrityksessä. Metsätieteen aikakauskirja 3/2002: 441–457. Packalén, P. & Maltamo, M. 2007. The k-MSN method for the prediction of species-specific stand attributes using airborne laser scanning and aerial photographs. Remote sensing of Environment 109(3): 328–341. Petterson, H. 1955. Barrskogens volymproduction. Med- delanden från Statens Skogsforskningsinstitut 45(1). 399 s. Prodan, M. 1965. Holzmesslehre. J.D. Sauerländers Ver- lag, Frankfurt Am Main. 644 s. Schmidt, M., Kiviste, A. & von Gadow, K. 2011. A spa- tially explicit height-diameter model for Scots pine in Estonia. European Journal of Forest Research 130: 303–315. doi10.1007/s10342-010-0434-8. Siipilehto, J. 1999. Improving the accuracy of predicted basal-area diameter distribution in advanced stands by determining stem number. Silva Fennica 33(4): 281–301. — 2000. A comparison of two parameter prediction meth- ods for stand structure in Finland. Silva Fennica 34(4): 331–349. — 2001. Improving the accuracy of predicted diameter and height distributions. Julkaisussa: LeMay, V. & Marshall, P. (toim.). Forest modelling for ecosystem management, forest certification, and sustainable management. Proceedings of the conference held in Vancouver, BC, Canada, August 12 to 17, 2001. s. 385–391. — 2009. Modelling stand structure in young Scots pine dominated stands. Forest Ecology and Management 257: 223–232. — 2011a. Local prediction of stand structure using linear prediction theory in Scots pine-dominated stands in 233 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… Finland. Silva Fennica 45(4): 669–692. — 2011b. Methods and applications for improving param- eter prediction models for stand structures in Finland. Dissertationes Forestales 124. http://www.metla.fi/ dissertationes/df124.htm — & Huttunen, T. 2015. Metsikön todellisen iän määrit- täminen rinnankorkeusiästä suomalaisissa talousmet- sissä. Ikälisäysmalli INKA-aineiston kairauksista. Metsätieteen aikakauskirja 2/2015: 87–104. — , Valkonen, S., Ojansuu, R., Hynynen J., Miina J. & Saksa, T. 2014. Metsikön varhaiskehityksen kuvaus MOTTI-ohjelmistossa. Metlan Työraportti 286. 43 s. Saatavissa: http://www.metla.fi/julkaisut/workingpa- pers/2014/mwp286.htm Sirkiä, S., Heinonen, J., Miina, J. & Eerikäinen, K. 2015. Subject-specific prediction using a nonlinear mixed model: consequences of different approaches. Forest Science 61(2): 205–212. Solmun maastotyöopas. 1997. Metsätalouden kehittämis- keskus Tapio. 82 s Staudhammer, C. & LeMay, V. 2000. Height prediction equations using diameter and stand density measures. Forestry Chronicle 76(2): 303–309. Suvanto, A., Maltamo, M., Packalén, P. & Kangas, J. 2005. Kuviokohtaisten puustotunnusten ennustaminen laserkeilauksella. Metsätieteen aikakauskirja 4/2005: 413–428. Uuttera, J, Anttila P., Suvanto A. & Maltamo M. 2002. Uudet kuvioittaisen arvioinnin menetelmät – Arvio soveltuvuudesta yksityismaiden metsäsuunnitteluun. (Tieteen tori). Metsätieteen aikakauskirja 3/2002: 523–531. — , Anttila P., Suvanto A. & Maltamo M. 2006. Yksityis- metsien metsävaratiedon keruuseen soveltuvilla kau- kokartoitusmenetelmillä estimoitujen puustotunnus- ten luotettavuus. Metsätieteen aikakauskirja 4/2006: 507–519. Temesgen, H. & Gadow, K.von. 2004. Generalized height-diameter models – an application for major tree species in complex stands of interior British Columbia. European Journal of Forest Research 123: 45–51. doi 10.1007/s10342-004-0020-z. Valtakunnan metsien 11. inventointi (VMI11). Maasto- työn ohjeet 2009. Metsäntutkimuslaitos. Saatavilla: http://www.metla.fi/ohjelma/vmi/vmi11-maasto-oh- je09-2p.pdf. Vastaranta, M., Ojansuu, R. & Holopainen, M. 2010. Puustotunnusten laskennallisen ajantasaistuksen luo- tettavuus – tapaustutkimus Pohjois-Savossa. Metsäti- eteen aikakauskirja 4/2010: 367–381. Veltheim, T. 1987. Puumallit männylle, kuuselle ja koi- vulle. Helsingin yliopisto. 60 s. Vestjordet, E. 1972. Diameterfordelinger og høydekur- ver for ensal- drede granbestand. Meddelelser fra det Norske Skogforsøksvesen 29: 482–557. Yuancai, L. & Parresol, B.R. 2001. Remarks on height- diameter modelling. USDA Forest Service. Southern Research Station. Research Note SRS-10. 5 s. Zhang, L. 1997. Cross-validation of non-linear growth functions for modelling tree height-diameter relation- ships. Annals of Botany 79: 251–257. 70 viitettä http://www.metla.fi/dissertationes/df124.htm http://www.metla.fi/dissertationes/df124.htm http://www.metla.fi/julkaisut/workingpapers/2014/mwp286.htm http://www.metla.fi/julkaisut/workingpapers/2014/mwp286.htm 234 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli Liite Esimerkissä esitetään pituuden sekamallin lokalisointi (kalibrointi) lineaarisen ennustamisen teo- riaan perustuen. Aineistoon sovitettu malli Näslundin pituuskäyrän parametreille oli muotoa: yij = xi0’b0 + b0i + (xi1’b1 + b1i)dij + eij jossa xi0 ja xi1 kuvaavat mallin kiinteän osan selittäjiä parametreille b0 ja b1 ja b0 ja b1 ovat niitä vastaavat kertoimet, b0i ja b1i ovat satunnainen metsikkö vakio ja satunnainen metsikkökerroin. Pituus ennustetaan uudistuskypsän männikön (INKA metsikkö 1194) tilavuustunnuksilla. Mallissa tarvittavat tunnukset olivat lämpösumma (DDY) 1176 °Cvrk, tukkipuun ja kuitupuun tilavuudet olivat (Tukki) 243 m3ha–1 (Kuitu) 39 m3ha–1 (ks. taulukko 7). Näillä tunnuksilla ennusteet ja vari- anssifunktion mukainen jäännösvarianssi männylle olivat: b0 = 3,128–0,537(DDY/1000) – 0,041ln(Tukki+2) – 0,414ln(Kuitu+2) = 0,733 b1 = 0,232 – 0,025ln(Tukki+2) + 0,023ln(Kuitu+2) = 0,180 s(e)2 = 0,2572(1000/DDY)0,914 = 0,057 Mallin satunnaisparametrien estimoidut keskihajonnat ja korrelaatiokerroin olivat: s(b0i) = 0,245, s(b1i) = 0,0220 ja corr(b0i, b1i) = –0,621. Metsikön lokalisointiin käytetyt viisi koepuuta olivat: Puut 1 2 3 4 5 Läpimitta, cm 19,4 21,1 24,0 23,8 24,4 Pituus, m 20,7 20,1 20,3 19,2 20,8 Ennustevirheet (y – Xa) lasketaan muunnoksen y ennusteesta ja mitatuista puista dj,hj: ŷ j = b̂0 + b̂1dj ja yj = dj hj −1,3( ) missä j = puut 1,..,5. Sekamalli voidaan esittää yleisessä muodossa: y = Xa + Zb + e Missä y on havaintovektori, X on matriisi selittävistä muuttujista, Z on mallimatriisi, a on kiinteiden parametrien vektori, b on tuntemattomien satunnaisparametrien vektori ja e on jäännösvirhe. Jos var(b) = D ja var(e) = R, niin y:n varianssikovarianssimatriisi voidaan esittää: var(y) = ZDZ’ + R 235 Siipilehto & Kangas Näslundin pituuskäyrä ja siihen perustuvia malleja läpimitan ja pituuden välisestä riippuvuudesta… Tuntemattomat satunnaisparametrit saadaan BLUP estimaatista (esim. Lappi 1991, Kangas ja Maltamo 2002): b = [Z’R–1Z + D–1]–1Z’R–1 (y – Xa) Missä (y – Xa) on vektori havaittujen ja ennustettujen y arvojen eroista. Kirjoitetaan edellä esitetyt matriisit kolmen mittaushavainnon esimerkille. Mallimatriisi Z satunnaiselle vakiolle ja satunnai- selle kertoimelle on: Z = 1 d1 1 d2 1 d3 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ja sen transpoosi, johon d:n arvot on sijoitettu ′Z = 1 1 1 19,4 21,1 24,0 ⎛ ⎝⎜ ⎞ ⎠⎟ Ennustevirheiden vektoriksi saatiin: y –Xa( ) = 0,181 0,337 0,455 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ Varianssit saadaan satunnaisparametrien estimoidun keskihajonnan neliönä, esim. var(b0i)=s(b0i)2 ja kovarianssi saadaan estimoidusta korrelaatiokertoimesta cov(b0i,b1i) = corr(b0i,b1i) × s(b0i) × s(b1i). Siten estimoitu varianssi-kovarianssimatriisi satunnaisparametreille on D: D = 0,060 −0,00335 −0,00335 0,000484 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ja sen käänteismatriisi on D−1 = 27,12 187,53 187,53 3363 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ Jäännösvirheen varianssi on diagonaalinen (var(e) = Rn×n) eli 0,055 R3×3. Sen käänteismatriisi saadaan varianssin 0,055 käänteisluvusta, joka on 18,18. Kun näistä elementeistä kirjoitetaan yhtälö b = [Z’R–1Z + D–1]–1Z’R–1 (y – Xa), saadaan metsikön satunnaisparametrit b0i ja b1i matriisista b: b = 1 19,4 1 21,1 1 24,0 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ′ 0,057 0 0 0 0,057 0 0 0 0,057 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ −1 1 19,4 1 21,1 1 24,0 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ + 0,060 −0,00335 −0,00335 0,000484 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ −1 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ −1 × 1 19,4 1 21,1 1 24,0 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ′ 0,057 0 0 0 0,057 0 0 0 0,057 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ −1 0,181 0,337 0,455 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ b = −0,047 0,016 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 236 Metsätieteen aikakauskirja 4/2015 Tutkimusartikkeli Ennustetut satunnaisparametrit sijoitetaan yhtälöön, joka on lokalisoidussa muodossa metsikölle i ja puulle j: yij = b0 + b0i + b1 dij + b1i dij = 0,733 – 0, 047 + 0,180 dij + 0,016dij = 0,686 + 0,196 dij ja puun pituus vastaavasti hij = ( dij /yij)2 + 1,3. 0 5 10 15 20 25 0 5 10 15 20 25 30 h, m d, cm INKA 1194 Kiinteä osa Lokalisoitu n = 1 Lokalisoitu n = 3 Lokalisoitu n = 5 Muut puut Liitekuva 1. Metsikön tilavuustunnuksilla ennustettu pituuskäyrä mallin kiinteän osan mukaan ja ensimmäisellä, kolmella ja viidellä mittaushavainnolla lokalisoituna. Puut tulivat mukaan mittausjärjestyksessä vaaleimmasta tummimpaan ympyrään. Metsikön muut mitatut koepuut on esitetty kolmioilla. _GoBack