Scientia Forestalis, volume 43, n. 107
p.627-637, setembro de 2015

Predição do volume de árvores integrando Lidar e Geoestatística

Predict volume of trees integrating Lidar and Geostatistics

Samuel de Pádua Chaves e Carvalho1
Luiz Carlos Estraviz Rodriguez2
Luciana Duque Silva3
Luis Marcelo Tavares de Carvalho4
Natalino Calegario4
Mariana Peres de Lima1
Carlos Alberto Silva5
Adriano Ribeiro de Mendonça6
Marcos Felipe Nicoletti7

1Doutor(a), Professor(a) Adjunto. UFMT - Universidade Federal de Mato Grosso - 78060-900 – Cuiabá, MT. E-mail: sam.padua@gmail.com; marianaperes@ufmt.br.
2Doutor, Livre Docente no Departamento de Ciências Florestais. USP – Universidade de São Paulo / ESALQ – Escola Superior de Agricultura “Luiz de Queiroz”. Av. Pádua Dias, nº 11 - Caixa Postal 09 – 13418-900 - Piracicaba, SP. E-mail: lcer@usp.br.
3Doutora, Professora no Departamento de Ciências Florestais. USP – Universidade de São Paulo / ESALQ – Escola Superior de Agricultura “Luiz de Queiroz”. Av. Pádua Dias, nº 11 - Caixa Postal 09 – 13418-900 - Piracicaba, SP. E-mail: luciana.duques@usp.br.
4Doutor, Professor Adjunto. UFLA - Universidade Federal de Lavras. Caixa Postal 3037- 37.200-000, Lavras, MG.  E-mail: passarinho@dcf.ufla.br; calegari@dcf.ufla.br.
5Mestre, Programa de Pós-Graduação em Recursos USP – Universidade de São Paulo / ESALQ – Escola Superior de Agricultura “Luiz de Queiroz”. Av. Pádua Dias, nº 11 - Caixa Postal 09 – 13418-900 - Piracicaba, SP. E-mail: carlos_engflorestal@outlook.com.
6Doutor, Professor Adjunto. UFES - Universidade Federal do Espirito Santo. Av. Governador Lindemberg, 316 - 29550-000, Jerônimo Monteiro-ES. E-mail: ribeiroflorestal@yahoo.com.br.
7Mestre, Professor do Departamento de Engenharia Florestal. UDESC – Universidade do Estado de Santa Catarina. Av. Luis de Camões, 2090 – 88.520-000 - Lages-SC. E-mail: nicoletti@cav.udesc.br.

Recebido em 10/05/2014 - Aceito para publicação em 23/03/2015

Resumo

Este estudo propõe integrar geoestatística, medições de circunferência em campo e escaneamento a laser para predição de volume de madeira. O método descrito nesse estudo considera a modelagem da variável circunferência das árvores por meio da geoestatística. Para este fim foram considerados dois cenários de ajustes: o primeiro considerou a krigagem ordinária e o segundo a cokrigagem ordinária, tendo como variável auxiliar, a altura das árvores, sendo estas oriundas dos sobrevoos LiDAR. De acordo com as estatísticas propostas, a cokrigagem ordinária foi significativamente superior à krigagem ordinária, em que o valor do Critério de Informação de Akaike (AIC) foi reduzido em 32,7 unidades, a raiz quadrada do erro médio (RMSE) em 40%, e 55% de superioridade para o coeficiente de determinação (r2), além da boa distribuição dos resíduos com médias próximas a zero. Após selecionado o modelo geoestatístico, foram geradas predições no gridde alturas gerado pelo processamento dos dados LiDAR. Após obtidos os pares de altura e circunferência, foi posteriormente aplicado o modelo de afilamento para gerar as predições dos volumes das árvores que compunham o talhão florestal. Os resultados permitiram concluir que o método proposto é tão preciso quanto aos levantamentos feitos por inventários florestais convencionais, com diferenças médias de 0,7% na estimativa do volume e 0,18% para número de árvores.
Palavras-chave: Continuidade Espacial; LiDAR; Inventário Florestal.

Abstract

This study aims to integrate spatial pattern modeled from field circumference measurements and airborne laser scanner data during volume estimation. The tree circumference determination was based in two approaches. In the first, the spatial variation of circumference is constant in average, and in the second, the spatial dependency of circumference was modeled based on the spatial distribution of height. The geo-statistical model considering spatial distribution of height was statistically superior based on Akaike's Information criterion, improving the performance in 32.7 units compared to the alternative modeling. Coefficient of determination also increased by 55%; no bias was detected, and the error was close to zero. The geo-statistical model estimated the circumference for trees extracted based on LiDAR data. Thus, the diameter and height was used as input to a logistic taper equation to estimate volume tree by tree. The results indicated that both methods showed similar results, differing by 0.7% as to volume and by 0.18% as to the number of trees.
Keywords: Spatial continuity; LiDAR; Forest Inventory.


INTRODUÇÃO

Em todos os ramos da ciência moderna é constante a busca por novas tecnologias capazes de fornecer informações com alto grau de precisão e rapidez. Na ciência florestal, em especial no Brasil, um exemplo recentemente empregado é o uso de tecnologias com equipamentos aero embarcados a laser, LiDAR (Light Detection and Ranging). Apesar da sua recente aplicação no setor florestal, seu uso já é corriqueiro em outros ramos da engenharia, como em projetos altimétricos e em projetos de redes viárias, ou ainda em transmissão de energia elétrica (RODRIGUEZ et al., 2010).

Na ciência florestal ressalta-se o pioneirismo do pesquisador Eric Naesset (NAESSET, 1997a, 1997b), com trabalhos originados inicialmente em países escandinavos. Segundo Nelson et al. (1984), é relevante a qualidade bem como o número de informações obtidas direta ou indiretamente nos levantamentos com o LiDAR. Pode-se dizer que informações úteis para aplicação em inventário florestal são subprodutos dos sobrevoos, uma vez que os principais produtos são o modelo digital do terreno e o modelo digital de superfície.

As principais informações utilizadas nos resultados de inventário florestal são o censo das alturas das árvores bem como a contagem de indivíduos. Essas duas, juntamente com os valores de circunferência ou diâmetro são algumas das principais informações a serem extraídas nos levantamentos LiDAR para quantificação e qualificação da biomassa florestal.

Mesmo que o censo da altura e a contagem dos indivíduos sejam as principais variáveis de obtenção nos sobrevoos a laser, outras importantes variáveis podem ser obtidas com coeficientes de correlação em torno de 70 e 80 % para área basal e biomassa aérea respectivamente (LEFSKY et al., 1999). Dean et al. (2009) associam a altura das árvores obtidas com o LiDAR com a altura de copa viva no dossel de florestas de Pinus no sudeste da Louisiana por meio da distribuição de probabilidade Weibull truncada. Os autores concluíram que, quando comparados às observações de campo, os erros não diferem.

Em se tratando do uso da ferramenta no Brasil, Macedo (2009) comparou dados de LiDAR aos mensurados em campo para plantios clonais de Eucalyptus na região do Vale do Paraíba – SP e encontrou nível de precisão superior a 90% para volume, altura, diâmetro e número de fustes e de 80% para área de copa. Investigações correlatas também foram abordadas em outros trabalhos como Pires (2005); Lingnau et al. (2007); Zandoná et al. (2008); Zonete et al. (2010); Giongo et al. (2010); Rodriguez et al. (2010); Oliveira et al. (2012); Silva et al. (2014), garantindo correlações, campo e LiDAR, superiores a 80%.

Conforme Hudak et al. (2002), poucos são os trabalhos que integram sensoriamento remoto a outras fontes de informações. Os autores tratam em seu trabalho a abordagem de integração de imagens de satélites e LiDAR para estimar altura de copa e citam que apesar do grande número de trabalhos publicados em ambos os temas, na sua grande maioria, aplicam as técnicas de maneira independente.

Outra abordagem de integração de técnicas é o do uso de métodos Geoestatísticos associados ao LiDAR. Tesfamichael et al. (2009), visando estimar o número de fustes por hectare por modelos de semivariograma com a finalidade de definir tamanho ótimo da janela para gerar o modelo de copa, concluiu que o tamanho ótimo da janela para plantios de Eucalyptus grandis no sul da África está entre os valores 2 e 5,4m. Segundo Finley et al. (2013), ignorar a dependência espacial, muitas vezes presente em observações de dados florestais, pode ocasionar incorretas inferências sobre os parâmetros do modelo, consequentemente sobre as predições geradas.

Diante do exposto, este trabalho parte da hipótese de que é possível integrar medições de campo, LiDAR e Geoestatística visando estimar o volume de madeira em nível de povoamento. Partindo-se desta hipótese, os objetivos deste trabalho foram: modelar a tendência espacial da variável circunferência considerando um modelo geoestatístico de média constante (krigagem ordinária) e outro que considera as variações da média como uma função da variável preditora altura total (cokrigagem ordinária); parametrização no processamento dos dados LiDAR visando a contagem de árvores e estimativa média da altura; e obtenção do volume total de madeira, todos em nível de povoamento.


MATERIAL E MÉTODOS


Área de estudo

O trabalho foi realizado em plantios clonais de Eucalyptus urograndis localizados no município de São Luiz do Paraitinga-SP.


Figura 1. Área de estudo.
Figure 1. Study area.

Quanto ao manejo empregado na área, todos os talhões estavam sob regime de alto fuste, em segundo ciclo, plantados unicamente com clones de Eucalyptus grandis e Eucalyptus urophylla, e espaçamento predominante de 3,0 m na entrelinha e 2,0 m na linha de plantio, o que totaliza 1667 covas por hectare.


Características e coleta dos dados

As medições em campo foram realizadas no mês de maio do ano de 2012. Neste instante os talhões se apresentavam com 5,4 anos. Ao todo foram mensuradas 1.745 árvores em 26 parcelas de formato circular e 400m² de área, distribuídas em três talhões que totalizavam 139 hectares de plantio.

A intensidade amostral adotada foi de uma parcela representando cinco hectares. Das parcelas foram mensurados: Circunferência à Altura do Peito (CAP), medida a 1,30 m do solo; e Altura Total (HT) de todas as árvores além da classificação qualitativa das mesmas. Foram ainda coletadas as coordenadas geográficas do centro das parcelas, com receptores geodésicos, garantindo precisão de aproximadamente 10 cm. Utilizou-se o sistema de projeção de Mercator (UTM) e sistema geodésico de referência SIRGAS 2000, sendo que toda a área de estudo encontra-se na zona 23 Sul.

Estatísticas descritivas de algumas das variáveis dendrométricas podem ser visualizadas na tabela 1.

Tabela 1. Estatísticas descritivas em nível de talhão para a área de estudo.
Table 1. Descriptive statistics in stand level.
CAP (cm) HT (m) HDOM (m) NFUSTES (ha-1) VOL (m³.ha-1)
Talhão 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
Mínimo 44,92 42,58 42,75 18,3 21,2 22,3 24,7 25,5 25,4 1469 1443 1349 286,1 237,66 258,54
Máximo 51,41 47,93 50,54 26,7 25,1 27,1 30,2 28,1 29,1 1700 1875 1725 390,25 361,06 377,75
Médio 46,89 46,06 47,24 23,8 23,8 23,9 27,3 26,3 26,7 1621 1627 1524 324,02 306,37 305,12
Desvio padrão 2,05 2,00 2,03 2,49 1,45 1,51 2,00 0,96 1,30 72,7 142,7 128,3 37,5 49,0 23,0
Coeficiente de Variação 4,39 4,36 4,31 10,4 6,11 6,36 7,36 3,68 4,90 4,49 8,77 8,42 11,5 16,0 7,57
CAP = circunferência a altura do peito; HT altura total; HDOM = altura média dominante segundo o conceito de Assman; NFUSTES = número de fustes por hectare; VOL = média do volume total com casca por hectare.

Os dados LiDAR foram gerados por meio de um vôo realizado no início de 2012 utilizando uma aeronave bimotor Piper Seneca II acoplada a uma câmera HARRIER 68i, IMU/GPS Applanix 510 e sensor LiDAR Riegl LMS Q680I. A tabela 2 apresenta mais detalhes relativos ao vôo.

Tabela 2. Detalhes do sobrevoo.
Table 2. Flightdetails.
Atributos Valores
Densidade de pulso laser 5 pontos/m²
Resolução espacial 0.5m
Bandas RGB e IR
Resolução espectral 32 bits
Velocidade 55 m/s (198 km/h)
Altura de voo 629,24 m
Ângulo de varredura 60º
Precisão 10 a 15 cm (Precisão da IMU)
Largura da faixa de cobertura 726,58 m
Taxa de leitura 200.000 pulsos
Frequência de varredura 300 KHz
Câmera HARRIER 68i
IMU/GPS Applainix 510
Datum SIRGAS 2000 / UTM Zona 23S


Softwares, processamento e parametrizações

O passo inicial do processamento dos dados do sobrevoo LiDAR consiste em sumarizar a nuvem de pontos para obtenção dos dados de altura. A altura foi gerada pela diferença entre o modelo digital de terreno e o modelo digital de superfície. Como produto final foram geradas informações tridimensionais sendo duas dimensões representadas pelos pares de coordenadas geográficas e a terceira dimensão pelos valores das alturas das árvores individualizadas. Todo esse procedimento de individualização das árvores foi feito no software FUSION versão 3.30. O diagrama da figura 2 tem por finalidade resumir os passos necessários para geração da nuvem de pontos sumarizada.


Figura 2. Processamento dos dados LiDAR para obtenção da altura.
Figure 2. Processing data LiDAR to determine the height.

O passo 1 consistiu em gerar um resumo detalhado da nuvem de pontos. O passo 2 transformou o Modelo Digital de Terreno (DTM) originalmente em formato raster com extensão “.asc” para o formato “.dtm” (formato especifico do FUSION) com 1 metro de resolução de pixel.. O passo 3 consistiu na normalização e recorte da nuvem de pontos com o shape dos talhões das áreas de interesse. Por fim, o quarto passo computou as métricas da nuvem de pontos LiDAR. Mais detalhes sobre as sintaxes das funções apresentadas neste trabalho podem ser encontrados em McGaughey (2010).

O procedimento de individualização das árvores que posteriormente formou o grid de predição dependeu fortemente das parametrizações do quarto passo do diagrama (Figura 2). Para as típicas florestas clonais de eucalipto plantadas no Brasil é recomendável que o tamanho do pixel seja proporcional ao espaçamento do plantio. Este valor deve estar próximo a área útil ocupada por cada árvore, sendo que para este estudo especificamente, este valor é de 6m². Recomendações semelhantes foram feitas por Tesfamichael et al. (2009). O grid de predição com as árvores individualizadas e com seus respectivos pares de coordenadas geográficas e altura é o produto das etapas descritas até este ponto.

Para compor as amostras que fariam parte do processo de ajuste dos modelos geoestatísticos foi aleatorizada uma árvore por parcela totalizando, portanto, 26 árvores com informações de coordenadas X e Y, CAP mensurada no campo e HT originada do sobrevoo. Optou-se pelo sorteio de uma árvore uma vez que a homogeneidade da floresta é uma característica marcante para os plantios monoclonais de eucalipto tipicamente plantados no Brasil, em que partiu-se do princípio que há uma alta intensidade amostral na mensuração de informações para as florestas extremamente homogêneas. Homogeneidade esta verificada na tabela 1, com os baixos valores dos coeficientes de variação para as diversas variáveis dendrométricas, como altura e diâmetro por exemplo.

Obtido os pares de circunferência e altura, foi posteriormente obtido o volume de cada árvore quem compuseram o grid de predição. As estimativas de volume foram geradas pela aplicação da função de afilamento, descrita pelo modelo 1 e proposta por Carvalho (2013).

  (1)

Em que: Vj é o volume da j-ésima árvore em m³; R varia de 0 ao rap, sendo rapj raio mensurado a altura do peito da j-ésima árvore em cm; htj é a altura obtida do processamento dos dados LiDAR da j-ésima árvore em m; é a constante pi igual a 3,14159265; Φi são os parâmetros da regressão.

O volume de cada talhão foi obtido pela soma dos valores dos volumes individuais de cada árvore. Mais detalhes sobre o uso e aplicação do modelo citado podem ser obtidos em Calegario (2002); Carvalho (2013).

Para fins de comparação e validação do método proposto neste estudo, as informações de campo foram processadas por meio das formulações da Amostragem Casual Simples. Todas as análises estatísticas foram feitas pelo software R, versão 2.15.3. Os principais pacotes utilizados nesse estudo foram: geoR; sp; spdep; maptools; gpclib; rgdal; proj4; raster e lattice, citados também por Diggle e Ribeiro Jr. (2007); Bisvand et al. (2008).


Ajuste e seleção dos modelos

Duas abordagens de ajuste foram consideradas para modelar as variações espaciais do CAP, assumindo as premissas básicas de que o processo a ser modelado é estocástico, estacionário e isotrópico, ou seja, é probabilístico com distribuição normal multivariada ou suas transformações por Box e Cox (1964) em que a variância é a mesma em todas as direções na área, independente da localização onde se mensurou o atributo, no caso o CAP.

A primeira abordagem assumiu as variações espaciais da variável CAP como sendo função apenas de sua posição geográfica, ou seja, ajusta-se um modelo geoestatístico com média constante, ou seja, é considerada apenas a estrutura de correlação existente entre as observações da variável analisada, CAP (krigagem ordinária).

A segunda abordagem assumiu que as variações espaciais da variável CAP é função não apenas de sua posição geográfica, mas também da sua altura total (cokrigagem). Este processo pode ser generalizado pela expressão 2.

(2)

Em que: S(x) é um processo estocástico gaussiano com valores de coordenadas x pertencentes ao conjunto dos números reais positivos. Todo processo dessa natureza é especificado segundo Diggle e Ribeiro Jr (2007) por uma função de média , por uma função de covariância  e por uma função de correlação , sendo u = x - x’.

Para as duas abordagens tratadas nesse estudo foram avaliados três modelos para ajustar os semivariogramas dos modelos exponencial, gaussiano e esférico, descritos pelas equações 5, 6 e 7 respectivamente. As funções de semivariância citadas por Hudak et al. (2002) para o modelo de média constante (krigagem ordinária) bem como para o modelo com média variável (cokrigagem ordinária) são descritas respectivamente pelas equações 3 e 4.

 (3)

Em que: γ(h) = semivariância em função do lag da distância h; N(h) = número de par de dados separados pelo lag da distância h; z = valor do atributo na posição uα e uα+h

   (4)

Em que: γij(h) = semivariância cruzada entre as variáveis i e j; zi = valor do i-ésimo atributo na posição uα e uα+h; zj = valor do j-ésimo atributo na posição uα e uα+h

 (5)

(6)

 (7)

Em que:Φ = parâmetro básico da função de correlação

Os três parâmetros dos semivariogramas são: nugget, sill e range. Nugget é expresso por  e denominado variância do nugget ou efeito pepita. Se >0 implica em processo espacial descontínuo. Sill corresponde à variância do processo que estamos modelando, ou mesmo a assíntota do semivariograma empírico denotado por . Range ou alcance prático, por definição e convenções geoestatísticas de acordo com Diggle e Ribeiro Jr. (2007), é a distância u0 em que , portanto  sendo que o variograma teórico de um processo espacial estocástico pode ser representado por  .

Proposto por Krige (1951), o método de predizer pontos não amostrados para modelos de média constante é a krigagem ordinária eq.(8). Krigagem é considerada como um procedimento de predição que minimize os erros sobre a média com pesos associados à distância dos pontos vizinhos, e para média em função de variáveis preditoras, cokrigagem ordinária eq. (9) a qual é uma extensão multivariada da krigagem ordinária que trata, não somente da autocorrelação espacial na variável de interesse, mas também da correlação cruzada entre a variável de interesse e a variável preditora (HUDAK et al., 2002).

 (8)

Em que: Z = variável primária ou de interesse; λα = peso; uα = localização.

(9)

Em que: λα1 = peso da variável primária ou de interesse; uα1 = localização da variável de interesse; λα2 = peso da variável secundária ou preditora; uα2 = localização da variável preditora.

Z segundo Diggle e Ribeiro Jr. (2000) apud Mello et al. (2006) tem distribuição normal multivariada de média μ e variância eq.(19).

 (10)

Onde μ é o vetor de média, no qual é constante na krigagem ordinária e  na cokrigagem ordinária para este estudo, e HT é altura total obtida do sobrevoo LiDAR. 1 é o vetor de valores 1; σ² é a variância populacional; é a matriz de correlação; I é a matriz identidade e  é a variância do erro aleatório.

A função de densidade de probabilidade de Z é definida pela expressão a seguir:

(11)

Em que:  : vetor da variável de interesse;  : vetor de média (constante na krigagem ordinária e função de variáveis preditoras – HT na cokrigagem ordinária);  : refere-se ao termo de variância do modelo multivariado, expresso por .

Como critérios de diagnóstico dos modelos foram calculadas as seguintes estatísticas: Critério de Informação de Akaike (AIC) eq.(12); Raiz Quadrada do Erro Médio (RMSE) eq.(13); Viés ou BIAS eq.(14); Coeficiente de Determinação eq.(15), além das análises gráficas de resíduos originadas da validação cruzada.


Critério de Informação de Akaike (AIC)

O Critério de Informação de Akaike estima a distância de Kullback-Leiber esperada entre dois modelos probabilísticos (AKAIKE, 1973). Esta distância relativiza modelos obtidos de amostras aleatórias à sua verdadeira distribuição de probabilidade.

(12)

Quanto menor o valor de AIC melhor o modelo ajustado. OAIC é amplamente utilizado como critério de seleção de modelo, uma vez que penaliza modelos com número excessivos de parâmetros, selecionando, portanto, modelos mais parcimoniosos. Não necessariamente os modelos avaliados deverão ser hierarquizados, basta apenas serem considerados como concorrentes.


Raiz quadrada do erro médio (RMSE)

Valor que representa a estimativa do desvio padrão e variância amostral, dado por:

(13)

Em que: Y = valor observado;  = valor predito; n = número de observações; p = número de parâmetros.

Quanto menor o valor da RMSE, melhor o modelo ajustado.


Viés ou Bias (B)

Valor que representa a média dos resíduos. Quanto mais próximo a zero menos tendencioso e preferível o modelo ajustado.

 (14)

Em que: B = Viés; Yi = CAP observado (cm); = CAP estimado (cm); N = número de observações.


Coeficiente de Determinação (R²)

De acordo com Schneider et al. (2009), o coeficiente de determinação permite medir o grau de explicação do modelo, medindo a proporção total de variação para a média.

  (15)

Em que: R² = coeficiente de determinação (%); = média amostral do CAP (cm).


RESULTADOS E DISCUSSÃO


Resultados do processamento dos dados de campo

Tabela 3. Estatísticas de posição e dispersão para o processamento da ACS da variável volume.
Table 3. Position and dispersion statistics to volume in Simple Casual Sample.
Talhão Número de Parcelas Média [m³.ha-1] Desvio Padrão [m³.ha-1] Coeficiente de Variação [%] Erro do Inventário [%] IC_Inferior [m³.ha-1] IC_Superior [m³.ha-1] IC_Inferior [m³] IC_Superior [m³]
1 9 324,0255 37,5464 11,6 8,9 295,2782 352,7727 13.550,31 16.188,73
2 6 306,3744 49,0283 16,0 16,7 255,0754 357,6734 10.307,59 14.453,58
3 11 305,1211 34,0926 11,2 7,5 282,3133 327,9288 14.866,61 17.268,73

Os resultados da tabela 3 indicam o alto grau de homogeneidade dos plantios estudados, com valores de coeficiente de variação do volume próximos a 11%, com exceção do talhão 2 que, por sua vez, é menos homogêneo. Esses valores são considerados baixos quando comparados, por exemplo, a plantios seminais de Pinus taeda no Brasil, que apresentam valores próximos a 25% (FERRAZ FILHO, 2009). Os resultados apresentados na tabela 3 foram as bases das discussões, comparações e validações do método proposto neste estudo.


Modelagem Geoestatística

De acordo com os resultados da tabela 4, o modelo Esférico para média constante é estatisticamente superior aos demais e, portanto mais indicado para a predição por krigagem ordinária. Entretanto para o modelo de média variável, o Exponencial foi superior aos demais. As curvas dos modelos propostos nos seus respectivos variogramas estão representadas nas figuras 3 e 4.

Tabela 4. Estatísticas de seleção dos modelosgeoestatísticos ajustados no estudo.
Table 4. Statistics of the selection for the fitted geostatistics models in this study.
Modelo Método Predição AIC RMSE (cm) BIAS (cm) R² (%)
Exponencial Krigagem Ordinária 192,5846 8,21 -0,2694 4,61
Esférico Krigagem Ordinária 191,5459 7,81 -0,3589 13,77
Gaussiano Krigagem Ordinária 191,7174 7,84 -0,3743 13,07
Exponencial CoKrigagem Ordinária 158,8363 4,64 -0,0062 69,44
Esférico CoKrigagem Ordinária 158,6754 4,65 0,0192 69,41
Gaussiano CoKrigagem Ordinária 158,3969 4,66 0,0417 69,27


Figura 3. Variograma empírico da variável CAP para o modelo de média constante (krigagem ordinária).
Figure 3. Empirical variogram of CAP variable to constant mean model.


Figura 4. Variograma empírico da variável CAP para o modelo de média variável (cokrigagem ordinária).
Figure 4. Empirical variogram of CAP variable to variable mean model.

Ao analisar separadamente (Tabela 5) os modelos considerados como superiores foi possível verificar a superioridade do modelo com média variável, utilizado na cokrigagem ordinária.

Tabela 5. Estatísticas de comparação do modelo Esférico para média constante versus o modelo Exponencial para média variável.
Table 5. Statistical comparison of the spherical model with constant mean versus exponential model with variavel mean.
Modelo Método de Predição AIC RMSE (cm) BIAS (cm) R² (%)
Esférico Krigagem Ordinária 191,5459 7,81 -0,3589 13,77
Exponencial CoKrigagem Ordinária 158,8363 4,64 -0,0062 69,44

As figuras 5 e 6 complementam os resultados encontrados anteriormente. A superioridade do modelo Exponencial é novamente comprovada em que se verifica uma maior concentração dos pontos na reta 45º conforme proposto por Graybill (1976).


Figura 5. Volumes preditos versus observados do modelo de média constante.
Figure 5. Predicted versus observed volume from constant mean model.


Figura 6. Volumes preditos versus observados do modelo de média variável.
Figure 6. Predicted versus observed volumes from variable mean model.

Os parâmetros dos modelos estão apresentados na tabela 6 nas escalas das transformações de Box e Cox (1964) sugeridas pelos valores de lambda utilizados para correção da normalidade, sempre que necessário como um procedimento padrão já implementado no pacote da geoR.

Tabela 6. Estimativa dos parâmetros de efeito da média (β), pepita ou nugget (), variância (),patamar (), alcance (), lambda (λ) e grau de dependência espacial (D.E.) dos modelos Esférico e Exponencial.
Table 6. Parameters estimate of mean (β), nugget (), variance(), patamar(), range(), lambda (λ), and degree of spatial dependence (D.E.) to spherical and exponential models.
Modelo D.E.(%)
Esferico 935,3403 24166 63041 87207 0,5448 1,91 72
Exponencial 2,6191 + 0,2248 0 0,106 0,106 0,0791 0,34 100


Krigagem e parâmetros dendrométricos.

O grid tridimensional gerado no processamento dos dados LiDAR, originado da nuvem de pontos sumarizada para extração da altura individual das árvores pode ser visualizado na figura 7.


Figura 7. Grid de predição para a variável altura total.
Figure 7. Grid of prediction for the total height variable.

Dentre as diversas métricas analisadas, a mais correlacionada com a altura das árvores foi o percentil 70 e, portanto utilizada para representar as alturas das árvores dos talhões em estudo. Zonete et al. (2010) ao modelarem atributos dendrométricos em plantios clonais de Eucalyptus no sul da Bahia, correlacionam a altura média por parcela com o percentil 10 e 90, com valores de R² superiores a 90%. Subentende-se que essas diferenças são devidas principalmente a dois aspectos: i) abordagens distintas, uma vez que os autores trabalham em nível de parcela; ii) a qualidade do isolamento das árvores está diretamente relacionada ao espaçamento de plantio, e neste caso tratamos esta consideração única e exclusivamente como um suposto teórico.

Os mapas de krigagem podem ser visualizados pelas figuras 8 e 9. É possível notar uma maior amplitude de variação para o mapa da cokrigagem ordinária. Esses resultados são coerentes com àqueles apresentados pela figura 7, grid de predição, uma vez que foi verificada heterogeneidade da variável preditora, altura total, que apesar de ser uma das mais homogêneas variáveis em plantios clonais de Eucalyptus, pode apresentar variabilidades muitas vezes não captadas pelas amostras de inventário em campo.


Figura 8. Mapa de krigagem ordinária para o modelo de média constante.
Figure 8. Ordinary kriging map for the constant mean model.


Figura 9. Mapa de cokrigagem ordinária para o modelo de média variável.
Figure 9. Ordinary cokriging map for the variable mean model.

Depois de gerada a krigagem, aplicou-se o modelo de afilamento (Eq.1) para predizer o volume individual das árvores que compuseram o grid de predição. O volume de cada talhão é representado pelo somatório dos volumes individuais das árvores deste grid. Os resultados comparativos das abordagens trazidas nesse estudo estão expressos na tabela 7.

Tabela 7. Resultado comparativo dos métodos de inventário, krigagem ordinária e cokrigagem ordinária na predição do volume total de madeira em nível de talhão.
Table 7. Comparative analysis for the forest inventory, ordinary kriging and cokriging for the volume prediction at stand level.
Método Talhão Área (ha) DAP (cm) HT (m) NFUSTES / ha VOL (m³ / ha) VT (m³) DIF VOL (%)
Inventário¹ 1 45,89 14,9 23,9 1621 324,03 14869,53 ---
Krigagem² 1 45,89 16,7 23,9 1586 372,47 17092,57 15,0
CoKrigagem² 1 45,89 15,3 23,9 1586 325,69 14946,11 0,5
Inventário¹ 2 40,41 14,6 23,9 1627 306,37 12380,59 ---
Krigagem² 2 40,41 16,7 23,3 1587 364,22 14717,94 18,9
CoKrigagem² 2 40,41 15,0 23,3 1587 307,79 12437,85 0,5
Inventário¹ 3 52,66 15,1 24,0 1524 305,12 16067,68 ---
Krigagem² 3 52,66 15,4 23,6 1593 318,48 16771,40 4,4
CoKrigagem² 3 52,66 14,7 23,6 1593 296,77 15627,90 -2,7
¹ - Resultados obtido por mensurações de campo
² - Resultados gerados pela nuvem de pontos LiDAR associada a medições do CAP
DAP – diâmetro medido a 1,30m do solo (cm); HT – altura total (m); NFUSTES / ha – número de fustes por hectare; VOL – volume total em metros cúbicos por hectare; VT – volume total do povoamento em metros cúbicos por hectare; DIF VOL – diferença percentual do volume gerado pelos diferentes métodos em relação ao volume do inventário obtido pela ACS.

Ao analisar os resultados da tabela 7 é notável novamente a superioridade do método de cokrigagem ordinária quando comparado ao método da krigagem ordinária, com diferenças superiores a 10%. O método da cokrigagem permitiu predizer valores muitos próximos quando comparado aos resultados dos dados de campo, com diferenças de 0,5 a 2,7%.

As diferenças médias ponderadas do volume foram de -0,7% para cokrigagem ordinária e 12,1% para krigagem ordinária. Os bons resultados obtidos na individualização das árvores estão diretamente associados às parametrizações indicadas no quarto passo do diagrama (Figura 2) e enfatizados por Tesfamichael et al., (2009), com diferenças médias de 0,18%. Ressalta-se ainda que toda análise comparativa teve como referência os valores médios do inventário sem levar em consideração a possibilidade de flutuação tanto para mais como para menos o que certamente possibilitaria inferir sobre erros ainda menores.

Visando eliminar o efeito da casualização, todos os procedimentos apresentados por essa abordagem foram simulados em mais nove amostras. Os resultados desta simulação estão demonstrados na tabela 8.

Tabela 8. Repetição do procedimento para nove amostras aleatórias.
Table 8. Procedure replication to nine random sample.
Amostra Vol_Inv (m³.ha-1) Vol_CoKrigagem (m³.ha-1) Diferença (%)
1 311,73 294,58 -5,5
2 311,73 313,76 0,7
3 311,73 305,51 -2,0
4 311,73 301,38 -3,3
5 311,73 305,69 -1,9
6 311,73 303,86 -2,5
7 311,73 303,60 -2,6
8 311,73 303,99 -2,5
9 311,73 298,34 -4,3
Vol_Inv = volume total de madeira gerado pela ACS; Vol_CoKrigagem = volume total de madeira gerado pela coKrigagem; Diferença = Vol_CoKrigagem / Vol_Inv.

Os resultados eliminam o efeito da casualização da amostra e comprova a eficiência do método proposto com diferenças médias inferiores a 6%.


CONCLUSÕES

Os resultados permitiram concluir que as hipóteses levantadas foram comprovadas, ou seja, o método é válido e permite integrar técnicas de sensoriamento remoto, LiDAR, à mensurações de campo, garantindo precisões muito próximas aos resultados que foram validados pelas mensurações de campo.

O modelo geoestatístico com média variável (cokrigagem ordinária) foi significativamente superior ao modelo de média constante (krigagem ordinária) sendo, portanto o mais indicado para explicar as variações espaciais da variável CAP e por consequência as estimativas volumétricas.

A baixa variabilidade das florestas clonais de Eucalyptus associada aos altos padrões de uniformidade destes plantios garantiu a boa acurácia nos resultados do isolamento das árvores, bem na geração das informações de altura total e na predição do CAP.


AGRADECIMENTOS

À Fibria Celulose pela concessão dos dados.


REFERÊNCIAS BIBLIOGRÁFICAS

AKAIKE, H. Information theory and an extension of the maximum likelihood principle. In: PETROV, B. N.; CSAKE, F. (Eds.) INTERNATIONAL SYMPOSIUM ON INFORMATION THEORY, 2, 1973, Budapest. Proceedings… Budapest: Akademiai Kiado, 1973. p. 267–281.

BISVAND, R. S.; PEBESMA, E. J.; GÓMEZ-RUBIO, V. Applied spatial data analysis with R. New York: Springer, 2008. 379 p.

BOX, G. E. P.; COX, D. R. An Analysis of Transformations. Journal of the Royal Statistical Society. Series B (Methodological), London, v. 26, n. 2, p. 211-252, 1964.

CALEGARIO, N. Modeling Eucalyptus stand growth based on linear and nonlinear mixed-effects models. 2002. 123 p. Tese (Doutorado em Ciências Florestais) – Universidade da Georgia, Athens, 2002. 

CARVALHO, S. P. C. Estimativa volumétrica por modelos mistos e tecnologia laser aerotransportado em plantios clonais de Eucalyptus sp. 2013. 104 p. Tese (Doutorado em Recursos Florestais) - Escola Superior de Agricultura “Luiz de Queiroz”, Universidade de São Paulo, Piracicaba, 2013.

DEAN, T. J.; CAO, Q. V.; ROBERTS, S. D.; EVANS, D. L. Measuring heights to crown base and crown median with lidar in a mature, even-aged loblolly pine stand. Forest Ecology and Management, Amsterdam, v. 257, n. 1, p. 126-133, 2009.

DIGGLE, P. J.; RIBEIRO JR., P. J. Model-based geostatistics. New York: Springer, 2007. 230 p.

FERRAZ FILHO, A. C. Sistema de prognose do crescimento e produção para Pinus taeda L. sujeito a regimes de debastes e podas. 2009. 158 p. Dissertação (Mestrado em Engenharia Florestal) – Universidade Federal de Lavras, Lavras, 2009.

FINLEY, A. O.; BANERJEE, S.; COOK, B. D.; BRADFORD, J. B. Hierarchical Bayesian spatial models for predicting multiple forest variables using waveform LiDAR, hyperspectral imagery, and large inventory datasets. International Journal of Applied Earth Observation and Geoinformation, v. 22, p.147-160, 2013.

GIONGO, M.; KOEHLER, H. S.; MACHADO, S. A.; KIRCHNER, F. F.; MARCHETTI, M. LiDAR: Princípios e aplicações florestais. Pesquisa Florestal Brasileira, Colombo, v. 30, n. 63, p. 231-244, 2010.

GRAYBILL, F. A. Theory and application of the linear model. Belmont: Duxbury Press, 1976. 720 p.

HUDAK, A. T.; LEFSKY, M. A.; COHEN, W. B.; BERTERRETCHE, M. Integration of lidar and Landsat ETM + data for estimating and mapping forest canopy height. Remote Sensing of Environment, New York, v. 82, n. 2-3, p. 397-416, 2002.

KRIGE, D. G. A statistical approach to some basic mine valuation problems on the Witwatersrand, Journal of the Chemical, Metallurgical and Mining Society of South Africa, Johannesburg,  v. 52, p. 119–139, 1951.

LEFSKY, M. A.; HARDING, D.; COHEN, W. B.; PARKER, G.; SHUGART, H. H. Surface lidar remote sensing of basal area and biomass in deciduous forest of eastern. Remote Sensing of Environment, New York, v. 67, n. 1, p. 83-98, 1999.

LINGNAU, C.; NAKAJIMA, N. Y.; DAMAS, B.; SANTOS, D. S.; VINHAL, L. A. Obtenção de parâmetros florestais através de laser terrestre – Novas perspectivas. In: SIMPÓSIO BRASILEIRO DE SENSORIAMENTO REMOTO, 13, 2007, Florianópolis, Anais... São José dos Campos: INPE, 2007. p. 3661-3663.

MACEDO, R. C. Estimativa volumétrica de povoamento clonal de Eucalyptus sp através de laserscanner aerotransportado. 2009. 145 p. Dissertação (Mestrado em Sensoriamento Remoto) – Instituto Nacional de Pesquisas Espaciais. São José dos Campos, 2009.

MCGAUGHEY, R. J. FUSION/LDV: Software for LIDAR Data Analysis and Visualization. United States Department of Agriculture. Forest Service. Pacific Northwest Research Station, Seattle, 150 p. Disponível em: . Acesso em: 02 jul. 2010.

MELLO, J. M.; OLIVEIRA, M. S.; BATISTA, J. L. F.; JUSTINIANO JR., P. R.; KANEGAE JR., H. Uso do estimador geoestatístico para predição volumétrica por talhão. Floresta, Curitiba, v. 36, n. 2, p. 251-260, 2006.

NAESSET, E. Determination of mean tree height of Forest stands using airborne laser scanner data. ISPRS Journal of Photogrammetry and Remote Sensing, Amsterdam, v. 52, n. 2, p. 49-56, 1997a.

NAESSET, E. Estimating timber volume of forest stands using airborne laser scanner data. Remote Sensing of Environment, New York, v. 61, n. 3, p. 246-253, 1997b.

NELSON, R.; KRABILL, W.; MACLEAN, G. Determining forest canopy characteristics using airborne laser data. Remote Sensing of Environment, New York, v. 15, n. 3, p. 201-212, 1984.

OLIVEIRA, L. T.; CARVALHO, L. M. T.; FERREIRA, M. Z.; OLIVEIRA, T. C. A.; ACERBI JR., F.W. Application of lidar to forest inventory for tree count in stands of Eucalyptus sp. Cerne, Lavras, v. 18, n. 2, p. 175-184, 2012.

PIRES, J. M. Uso do LiDAR (Light Detection and Ranging) para estimação da altura de árvores em povoamentos de eucalipto. 2005.  50 p. Dissertação (Mestrado em Ciências Florestais) – Universidade Federal de Viçosa. Viçosa, 2005. 

RODRIGUEZ, L. E. C.; POLIZEL, J. L.; FERRAZ, S. B.; ZONETE, M. F.; FERREIRA, M. Z. Inventário florestal com tecnologia laser aerotransportada de plantios de Eucalyptus spp no Brasil. Ambiência, Guarapuava, v. 6, Ed. Especial, p. 67-80, 2010.

SCHNEIDER, P. R.; SCHNEIDER, P. S. P.; SOUZA, C. A. Análise de regressão aplicada à engenharia florestal. Santa Maria: Facos, 2009. 294 p.      

SILVA, C. A; KLAUBERG, C.; CARVALHO, S. P. C.; HUDAK, A. T.; RODRIGUEZ, L. C. E. Mapping aboveground carbon stocks using LiDAR data in Eucalyptus spp. plantations in the state of São Paulo, Brazil. Scientia Forestalis, Piracicaba, v. 42, n. 104, p. 200-215, dez. 2014.

TESFAMICHAEL, S. G.; AHMED, F.; VAN AARDT, J. A. N.; BLAKEWAY, F. A semi-variogram approach for estimating stems per hectare in Eucalyptus grandis plantations using discrete-return lidar height data. Forest Ecology and Management, Amsterdam, v. 258, n. 7, p. 1188-1199, 2009.

ZANDONÁ, D. F.; LINGNAU, C.; NAKAJIMA, N. Y. Varredura a laser aerotransportado para estimativa de variáveis dendrométricas. ScientiaForestalis, Piracicaba, v. 36, n. 80, p.295-306, dez. 2008.

ZONETE, M. F.; RODRIGUEZ, L. C. E.; PACKÁLEN, P. Estimação de parâmetros biométricos de plantios clonais de eucalipto no sul da Bahia: uma aplicação da tecnologia laser aerotransportada. Scientia Forestalis, Piracicaba, v. 38, n. 86, p. 225-235, jun. 2010.