Scientia Forestalis, volume 43, n. 107
p.601-608, setembro de 2015

Propostas de mapeamentos da capacidade produtiva de sítios florestais por meio de análises geoestatísticas

Proposals to mapping of forest site productivity by geo-statistical analyzes

Allan Libanio Pelissari1
Sidney Fernando Caldeira2
Afonso Figueiredo Filho3
Sebastião do Amaral Machado4

1Bolsista CNPq - Brasil. Doutor em Engenharia Florestal. UFPR - Universidade Federal do Paraná. Av. Pref. Lothário Meissner, 900 Jardim Botânico - 80210170 - Curitiba, PR. E-mail: allanpelissari@gmail.com.
2Professor Associado do Departamento de Engenharia Florestal. UFMT -  Universidade Federal de Mato Grosso. Av. Fernando Correa da Costa s/n Coxipó 78068-000 - Cuiaba, MT. E-mail: sidneycal@gmail.com.
3Professor Sênior do Departamento de Engenharia Florestal. UNICENTRO - Universidade Estadual do Centro-Oeste - Centro de Ciências Agrárias e Ambientais. PR-153 - KM 7 – Riozinho - 84500-000 - Irati, PR. E-mail: afigfilho@gmail.com.
4Professor Sênior do Departamento de Ciências Florestais. UFPR - Universidade Federal do Paraná. Rua Lothário Meissner, 632 - Jardim Botânico 80210-170 - Curitiba, PR. E-mail: profsamachado@gmail.com.

Recebido em 09/07/2014 - Aceito para publicação em 03/03/2015

Resumo

O trabalho objetivou aplicar análises geoestatísticas para modelar e mapear a variabilidade espacial da capacidade produtividade de sítios florestais. A partir de parcelas permanentes alocadas em plantios de teca, foram determinadas as curvas de sítios de produtividades alta, média e baixa, tendo o décimo segundo ano com a idade de referência. Nesse período, a krigagem ordinária foi utilizada para modelar os padrões espaciais da altura dominante, considerando as classes de sítio nos mapas temáticos. Adicionalmente, com a krigagem indicatriz foram determinadas as probabilidades de ocorrência de sítios mais produtivos, por meio da transformação binária da altura dominante, com o sítio de qualidade média como valor de corte. Nos semivariogramas ajustados foi verificada a dispersão branda das observações em torno das estimativas, com o início crescente das semivariâncias em função da distância e a posterior estabilização. Além disso, a característica isotrópica foi corroborada pela semelhança estrutural dos semivariogramas direcionais. Por meio da krigagem ordinária, os limites espaciais das classes de sítio foram identificados e mapeados, ao passo que com a krigagem indicatriz foram delimitados os locais com a probabilidade maior de obter sítios mais produtivos.
Palavras-chave: Classes de produtividade; Variabilidade espacial; Krigagem indicatriz.

Abstract

The study aimed to apply geo-statistical analyses to model and to map the spatial variability of forests sites productivity. Through permanent plots allocated in teak stands, sites curves of high, medium and low yields were obtained, with the twelfth year as reference age. In this period, ordinary kriging was used to model the spatial patterns of the dominant height, considering the site classes in thematic maps. Additionally, through the indicator kriging highest sites probability was determined, with the binary transformation of dominant height, being the site of average quality cutoff value. In the adjusted semi-variogram a bland dispersion of observations around the estimates was verified, beginning with the increase of the semi-variance versus distance and the subsequent stabilization. In addition, the isotropic characteristic was confirmed by the structural similarity of the directional semi-variograms. Through ordinary kriging, the spatial limits of the site classes were identified and mapped, while, with the indicator kriging, locals were delimited, which had the highest probability of being highly productive sites.
Keywords: Productivity classes; Spatial variability; Indicator kriging.


INTRODUÇÃO

A aparente homogeneidade espacial da estrutura dos povoamentos florestais é um aspecto que, frequentemente, dificulta observar as significativas variações existentes ao longo das áreas florestadas (MELLO et al., 2005a; RUFINO et al., 2006). Com isso, por meio dos avanços tecnológicos da silvicultura de precisão é possível mensurar e identificar as relações espaciais dos fatores que limitam a produtividade dos plantios, principalmente com os métodos geoestatísticos aplicados na predição de valores em locais não amostrados.

Na Teoria das Variáveis Regionalizadas (MATHERON, 1971), a geoestatística fundamenta a variável regionalizada como uma função espacial numérica de um fenômeno caracterizado no espaço, em que a semivariância é a sua medida estatística básica, por meio da qual é mensurada a estrutura espacial e as relações estatísticas existentes entre pontos amostrais separados por sucessivas distâncias (CARVALHO; VIEIRA, 2001; DAVIS, 2002; ABREU et al., 2003).

Embora seja amplamente aplicada na geologia e na ciência do solo, a geoestatística possui potencial para descrever o comportamento espacial de variáveis dendrométricas de espécies florestais (NANOS et al., 2004; MELLO et al., 2005a; RUFINO et al., 2006; PEREIRA et al., 2011; ROSA FILHO et al., 2011; PELISSARI et al., 2012), algo não identificável pela estatística clássica e que, por conseguinte, acarreta em perda de informação (NANOS et al., 2005; AMARAL et al., 2010; LEAL et al., 2011).

Nesse contexto, como alternativa às análises espaciais tradicionais, a krigagem indicatriz é um método geoestatístico não paramétrico que utiliza a posição e os valores de dados codificados em zero ou um, em função de um valor de corte previamente estabelecido, visando o mapeamento da probabilidade de ocorrência das propriedades espaciais de um determinado fenômeno (MOTOMIYA et al., 2006; ASSUMPÇÃO et al., 2007; PAZ-FERREIRO et al., 2010).

Dessa forma, com o conhecimento técnico-científico das características espaciais dos maciços florestais é possível recomendar intervenções localizadas e tratos culturais e silviculturais direcionados para as diferentes condições locais. Dessa forma, este trabalho objetivou aplicar e avaliar análises geoestatísticas para a modelagem da variabilidade espacial e o mapeamento da capacidade produtividade do sítio em povoamentos florestais.


MATERIAL E MÉTODOS


Local de estudo e determinação das classes de sítio

O estudo foi desenvolvido em povoamentos de teca implantados em 1.260 hectares no espaçamento 3 m x 3 m e localizados no estado de Mato Grosso sob as coordenadas geográficas 16°13'30'' S a 16°13'50'' S e 56°22'30'' W a 56°24'30'' W. O clima da região foi classificado como Aw (Köppen), com precipitação média de 1.300 a 1.600 mm ano-1 e temperatura média anual de 24°C a 26°C (ALVARES et al., 2013), e o solo foi identificado como Planossolo Háplico eutrófico de textura franco-argilo-arenosa em relevo suavemente ondulado (EMBRAPA, 2006).

Foram utilizados dados provenientes de uma série histórica, com idades entre dois a doze anos (Tabela 1), de 273 parcelas permanentes georreferenciadas de 15 m x 30 m. Nessas unidades amostrais foram obtidos os valores da altura dominante () de Assmann (1970), a qual expressa as interações dos fatores bióticos, edáficos e climáticos que definem a capacidade produtividade de um sítio florestal, sendo, ainda, a variável dendrométrica não influenciada significativamente pelos tratos silviculturais.

Tabela 1. Análise estatística descritiva da altura dominante nos povoamentos de teca.
Table 1. Descriptive statistical analysis of the dominant height in teak stands.
Idade
(ano)
Mínimo Média Máximo Desvio padrão Coeficiente de variação Grubbs Kolmogorov
-Smirnov
(m)
2 3,13 5,57 7,88 0,99 17,81% 2,467NS 0,049NS
3 5,65 9,43 12,23 1,48 15,69% 2,556NS 0,083NS
4 6,73 11,31 13,45 1,18 10,41% 2,897NS 0,032NS
5 10,63 12,92 15,05 0,92 7,13% 2,491NS 0,031NS
6 10,90 14,52 17,28 1,01 6,92% 2,747NS 0,092NS
7 14,00 16,63 19,97 0,91 5,50% 2,881NS 0,094NS
8 15,43 18,31 21,55 1,08 5,88% 3,003NS 0,091NS
9 16,76 19,00 20,53 0,62 3,24% 2,641NS 0,085NS
10 17,22 20,00 21,85 0,70 3,50% 2,444NS 0,087NS
11 18,65 22,00 24,16 0,89 4,06% 2,448NS 0,084NS
12 20,73 22,52 24,19 0,75 3,35% 2,374NS 0,083NS
Para teste Grubbs: NS = não há valores outliers na série de dados; e
Para teste Kolmogorov-Smirnov (KS): NS = há distribuição normal.

Por meio do prévio ajuste do modelo não linear Chapman (1961) e Richards (1959), a equação  foi utilizada para estimar a altura dominante em função da idade, cujos índices de qualidade de ajustamento apresentaram coeficiente de determinação ajustado igual a 0,95 e erro padrão da estimativa em porcentagem de 7,67%, além da ausência de resíduos tendenciosos (Figura 1A).


Figura 1. Distribuição dos resíduos (A) e curvas de sítio (B) obtidas pelo modelo Chapman-Richards nos povoamentos de teca.
Figure 1. Residues distribution (A) and site curves (B) obtained by Chapman-Richards model in teak stands.

Com base nessa equação ajustada e no método da curva-guia (SCOLFORO, 2006), curvas de sítio foram construídas para a idade de referência (IREF) de 12 anos, uma vez que correspondeu ao período mais próximo à rotação técnica da cultura com observações coletadas. Com isso, três classes de sítio foram obtidas, cujos valores das alturas dominantes expressaram os índices de sítio. Dessa forma, com os sítios de alta (I), média (II) e baixa (III) produtividades foram obtidos os índices de sítios de 26, 22 e 18 m, respectivamente (Figura 1B).


Análises geoestatísticas da capacidade produtiva do sítio

A análise geoestatística foi utilizada para modelar os padrões espaciais da altura dominante na idade de referência de 12 anos, por meio do cálculo das semivariâncias (1), considerando o posicionamento geográfico central das unidades amostrais no campo e o posterior cômputo das distâncias (h) e das diferenças numéricas da variável na malha de pontos.

(1)

Em que:  = semivariância da variável ;  = distância; e  = número de pares de pontos medidos  e , separados por .

Adicionalmente, para a determinação da probabilidade de ocorrência de sítios mais produtivos ao longo da área florestada, procedeu-se a transformação dos dados de altura dominante em zero ou um na idade de referência de 12 anos dos povoamentos, sendo respectivamente acima ou abaixo do valor de corte de 22 m, estabelecido como o centro de classe do sítio de qualidade média. Com isso, a transformação resultou em um conjunto de dados binários, o qual foi submetido à análise da continuidade espacial por meio do semivariograma indicador (2):

(2)

Em que:  = semivariância da variável ;  = distância;  = número de pares de pontos medidos  e , separados por ; e  = valor de corte.

As semivariâncias foram determinadas entre os pontos amostrais equidistantes, com a regularização da malha amostral por meio de uma tolerância angular de 22,5°, passo de 300 m e largura máxima de 3.000 m (FIGURA 9B), de modo a obter semivariogramas com maior número de pares de dados e mais suavizados (YAMAMOTO; LANDIM, 2013). Esse processo foi repetido em quatro direções, 0° (S-N); 45° (SO-NE); 90° (L-O); e 135° (NO-SE), dos quais foram obtidas as semivariâncias médias entre distâncias equivalentes e verificado a ausência de anisotropia.

Para as estimativas das semivariâncias, os modelos esférico (3), exponencial (4) e gaussiano (5) foram testados, cuja inclinação da tangente à origem  no modelo esférico é igual a , a qual corta o patamar no ponto em que ; ao passo que o alcance no exponencial tem significado analítico, uma vez que o patamar é alcançado quando , porém utiliza-se na prática igual a ; contudo o alcance no modelo gaussiano é definido como sendo  e o efeito pepita () usualmente indica erros de mensuração dos dados (DAVIS, 2002; KANEVSKI; MAIGNAN, 2004; WEBSTER; OLIVER, 2007).

(3)

(4)

(5)

Em que:  = semivariância da variável ;  = distância;  = efeito pepita;  = variância a priori; e = alcance.

A estrutura do semivariograma teórico foi composta pelo efeito pepita (), correspondente ao valor da semivariância para a distância zero; o patamar (), que representa a estabilização dos valores do semivariograma próxima à variância dos dados; a variância a priori (), dada pela diferença entre o patamar () e o efeito pepita (); e o alcance (), definido pela distância onde o semivariograma alcança o patamar, indicando o limite onde as unidades amostrais são correlacionadas entre si (VIEIRA, 2000).

Por meio de planilhas eletrônicas com função de otimização, o método dos mínimos quadrados ponderados (MELLO et al., 2005b; AZEVEDO et al., 2012) foi utilizado para os ajustes dos semivariogramas teóricos, visando minimizar a soma de quadrados dos desvios ponderados (SQDP), onde as diferenças quadráticas entre as semivariâncias observadas e as estimadas foram ponderadas de acordo com o número de pares de pontos utilizados para o cálculo médio das semivariâncias observadas em cada distância que compõe o semivariograma.

A avaliação e a seleção dos melhores ajustes dos semivariogramas teóricos foram baseadas na menor soma de quadrados dos desvios ponderados (SQDP), no maior coeficiente de determinação (R2) e na validação cruzada, a qual, quando ideal, fornece o coeficiente linear igual a zero e os coeficientes angular e de determinação da validação cruzada (R2vc) iguais a um.

A interpolação e a espacialização da altura dominante foram realizadas por meio da krigagem, a qual considera a dependência espacial e estima sem tendência e com variância mínima para a confecção de mapas temáticos (ANDRIOTTI, 2003; YAMAMOTO; LANDIM, 2013), sendo esses elaborados com o programa SURFER 11.0 versão demonstração (GOLDEN SOFTWARE, 2002).

Mais precisamente, foi utilizada a krigagem ordinária pontual com a geração de uma grade virtual de pontos amostrais regularmente espaçados em 50 m. Dessa forma, com os parâmetros obtidos dos ajustes dos semivariogramas e com os valores observados das unidades amostrais vizinhas, a altura dominante foi estimada nos pontos não amostrados na área florestada, por meio da formulação (6), considerando os valores das alturas dominantes dos limites das classes de sítio nos mapas temáticos.

(6)

Em que:  = estimador de krigagem;  = peso;  = dados experimentais; e  = número de dados.

Ademais, por meio da krigagem indicatriz (7), foram gerados os mapas temáticos da probabilidade de ocorrência de sítios mais produtivos, cujas estimativas resultaram nos valores de incerteza espacial, entre zero e 100%, esperados nos locais com menores e maiores produtividades, respectivamente, considerando a altura dominante de 22 m da classe de média produtividade como o valor de corte.

(7)

Em que:  = indicador;  = números de vizinhos;  = peso; e  = dado transformado em um indicador.


RESULTADOS E DISCUSSÃO

Com os valores de altura dominante das unidades amostrais na idade de referência de 12 anos dos povoamentos de teca, foi observada variabilidade de 3,35% pelo coeficiente de variação, com média aritmética de 22,52 m, além de distribuição normal pelo teste Kolmogorov-Smirnov e ausência de valores discrepantes pelo teste Grubbs, ambos ao nível de 5% de probabilidade (Tabela 1).

Por meio do efeito pepita (), que representou a variância não identificada nos semivariogramas da altura dominante, foram verificados valores inferiores a uma unidade (Tabela 2), ao passo que os alcances () representaram as distâncias máximas em que dois pontos amostrais correlacionaram-se espacialmente (REICHERT et al., 2008), enquanto as determinações em distâncias superiores retrataram independências entre si (WEBSTER; OLIVER, 2007).

Tabela 2. Parâmetros dos semivariogramas ajustados para a altura dominante nos povoamentos de teca.
Table 2. Semi-variograms parameters adjusted for dominant height in teak stands.
Modelo C0 C α (m) R2 SQDP
Krigagem ordinária
1 Esférico 0,517 0,338 2.896 0,958 3,6 × 10-4
2 Exponencial 0,381 0,473 2.506 0,953 5,4 × 10-4
3 Gaussiano 0,572 0,287 2.518 0,925 4,9 × 10-4
Krigagem indicatriz
1 Esférico 0,189 0,052 1.200 0,850 2,9 × 10-5
2 Exponencial 0,157 0,084 1.000 0,878 2,7 × 10-5
3 Gaussiano 0,188 0,053 800 0,839 3,4 × 10-5
Em que: C0 = efeito pepita; C = variância a priori; α = alcance;
R² = coeficiente de determinação; e SQDP = soma de quadrados dos desvios ponderados.

Além disso, como os valores dos coeficientes de determinação (R2) foram superiores a 0,92 e as somas de quadrados dos desvios ponderados (SQDP) variaram entre 3,6 × 10-4 a 5,4 × 10-4 (Tabela 2) para a krigagem ordinária, foi constatada a eficácia do modelo esférico para descrever o padrão espacial da altura dominante, para a qual foram obtidos o valor menor de SQDP e o maior de R2.

Ademais, com os valores de altura dominante, procedeu-se a transformação binária dos dados, considerando o valor de corte igual a 22 m na idade de referência de 12 anos, o qual correspondeu ao centro de classe do sítio de produtividade média. Dessa forma, procedeu-se o ajuste dos semivariogramas indicativos e, com base na menor soma de quadrados dos desvios ponderados (SQDP) e no maior coeficiente de determinação (R2), o modelo exponencial foi o selecionado (Tabela 2).

Na validação cruzada dos ajustes elegidos para as modelagens espaciais da altura dominante (Tabela 3), foram observados coeficientes lineares e angulares mais próximos aos ideais apenas na krigagem indicatriz, ao passo que os coeficientes de determinação (R2vc) foram baixos e os erros padrões de estimativas (Syx%) inferiores a 10% para duas modalidades de krigagem. Todavia, com a validação cruzada foi observado que os ajustes não foram necessariamente incorretos, visto que esses resultados são comumente encontrados em modelagem espaciais (ANDRIOTTI, 2003).

Tabela 3. Parâmetros da validação cruzada dos semivariogramas selecionados para a altura dominante nos povoamentos de teca.
Table 3. Cross-validation parameters of semi-variograms selected for dominant height in teak stands.
Modelo selecionado Coeficiente vc Syx%
Linear Angular
Krigagem ordinária
1 Esférico 8,117 0,373 0,344 5,80
Krigagem indicatriz
2 Exponencial 0,092 0,792 0,379 9,47
Em que: R²vc = coeficiente de determinação da validação cruzada; e Syx% = erro padrão da estimativa em porcentagem.

Visualmente, com os semivariogramas selecionados foram verificados os espalhamentos brandos das observações em torno das estimativas (Figuras 2A e 2B), com o início crescente das semivariâncias estimadas em função da distância e a posterior estabilização. Além disso, a característica isotrópica das modelagens foi corroborada pela semelhança estrutural dos semivariogramas direcionais (Figuras 2C e 2D).


Figura 2. Semivariogramas teóricos (A e B) e direcionais (C e D) ajustados para a altura dominante nos povoamentos de teca.
Figure 2. Theoretical (A and B) and directional semi-variograms (C and D) adjusted for dominant height in teak stands.

Com isso, após a seleção dos ajustes (Tabela 1) e constatada a dependência espacial entre as unidades amostrais e a ausência de anisotropia (Figura 2), procedeu-se a interpolação da altura dominante e o mapeamento das classes de sítio (Figura 3A) e das probabilidades de sítios mais produtivos nos povoamentos de teca (Figura 3B), respectivamente mediante a krigagem ordinária pontual e a krigagem indicatriz.


Figura 3. Mapas temáticos da distribuição espacial das classes de sítio (A) e das probabilidades de sítios mais produtivos (B) nos povoamentos de teca.
Figure 3. Spatial distribution thematic maps of site classes (A) and high site productive probabilities (B) in teak stands.

Por meio do mapeamento do sítio (Figura 3A) foram observadas delimitações espaciais definidas das classes de produtividade, com áreas aproximadas de 501; 700; e 59 hectares para os sítios de classes I, II e III, respectivamente. Dessa forma, foi constatada a eficácia da modelagem geoestatística para o zoneamento da capacidade produtiva do local, o que possibilita buscar evidências de características do meio que possam restringir o desenvolvimento da teca, tal como as limitações das propriedades físico-químicas do solo para a espécie, a fim de indicar tratos culturais adequados para o manejo desses povoamentos.

Além disso, por meio dos valores das probabilidades (Figura 3B) foram delimitadas as áreas com a possibilidade de obter sítios mais produtivos, o que propicia, dessa maneira, planejar a estrutura e a condução dos plantios, uma vez que os desbastes poderão apresentar maiores frequências ou intensidades nas áreas de qualidade alta, em decorrência do desenvolvimento superior dos indivíduos, enquanto nos locais de produtividade inferior, a densidade inicial dos plantios poderá ser menor, de modo a assegurar a produção sustentada e a mitigação de impactos negativos no meio, principalmente sob as reservas minerais do solo.

Ainda, ao considerar o conceito de qualidade de um sítio florestal como a soma das interações dos fatores bióticos, edáficos e climáticos que limitam a produtividade de um local (SPURR, 1952; CLUTTER et al., 1983), o emprego das análises espaciais permitiu que os aspectos das associações desses fatores do meio e das relações entre as unidades amostrais não fossem ignoradas e, desse modo, houve precisão estatística na composição dos mapeamentos das classes de sítio.


CONCLUSÕES

Por meio da krigagem ordinária aplicada à altura dominante dos povoamentos de teca, os limites espaciais das classes de sítio são identificados e mapeados ao longo das áreas florestadas, ao passo que com a krigagem indicatriz são delimitados os locais com a probabilidade de obter sítios mais produtivos para o planejamento da estrutura e condução dos plantios.


AGRADECIMENTOS

À empresa Teca do Brasil Ltda., pela disponibilização do banco de dados e pelo apoio técnico e logístico.


REFERÊNCIAS BIBLIOGRÁFICAS

ABREU, S. L.; REICHERT, J. M.; SILVA, V. R.; REINERT, D. J.; BLUME, E. Variabilidade espacial de propriedades físico-hídricas do solo, da produtividade e da qualidade de grãos de trigo em Argissolo Franco Arenoso sob plantio direto. Ciência Rural, Santa Maria, v. 33, n. 2, p. 275–282, 2003.

ALVARES, C. A.; STAPE, J. L.; SENTELHAS, P. C.; GONÇALVES, J. L. M.; SPAROVEK, G. Köppen’s climate classi?cation map for Brazil. Meteorologische Zeitschrift, v. 22, p. 1-18, 2013.

AMARAL, L. P.; FERREIRA, R. A.; WATZLAWICK, L. F.; GENÚ, A. M. Análise da distribuição espacial de biomassa e carbono arbóreo acima do solo em floresta ombrófila mista. Ambiência, Guarapuava, v. 6 (Especial), p. 103–114, 2010.

ANDRIOTTI, J. L. S. Fundamentos de estatística e geoestatística. São Leopoldo: UNISINOS, 2003. 165 p.

ASSMANN, E. The principles of forest yield study: studies in the organic production, structure, increment, and yield of forest stands. Oxford: Pergamon Press, 1970. 506 p.

ASSUMPÇÃO, R. A. B.; URIBE-OPAZO, M. A.; SOUZA, E. G.; JOHANN, J. A. Uso da krigagem indicatriz na avaliação da probabilidade da produtividade de soja segundo os padrões regional, estadual e nacional. Acta Scientiarum. Agronomy, Maringá, v. 29, n. 2, p. 165–171, 2007.

AZEVEDO, C. A. V.; PORDEUS, R. V.; DANTAS NETO, J.; AZEVEDO, M. R. Q. A.; LIMA, V. L. A. Dependência espacial da qualidade da água subterrânea no perímetro irrigado de São Gonçalo, Paraíba. Revista Verde de Agroecologia e Desenvolvimento Sustentável, v. 7, n. 2, p. 129–136, 2012.

CARVALHO, J. R. P.; VIEIRA, S. R. Avaliação e comparação de estimadores de krigagem para variáveis agronômicas – uma proposta. Campinas: Embrapa Informática Agropecuária, 2001. 21 p.

CHAPMAN, D. G. Statistical problems in dynamics of exploited fisheries populations. In: NEYMAN, J. (Org.) BERKELEY SYMPOSIUM ON MATHEMATICAL STATISTICS AND PROBABILITY: CONTRIBUTIONS TO BIOLOGY AND PROBLEMS OF MEDICINE, 4, 1961, Berkeley, Proceedings… Berkeley: University of California Press, 1961, v. 4, p. 153–168.

CLUTTER, J. L.; FORTSON, J. C.; PIENAAR, L. V.; BRISTER, G. H.; BAILEY, R. L. Timber management: a quantitative aprroach. New York: John Wiley & Sons, 1983. 333 p.

DAVIS, J. C. Statistic and data analysis in geology. 3.ed. New York: John Wiley & Sons, 2002. 656 p.

EMBRAPA - EMPRESA BRASILEIRA DE PESQUISA AGROPECUÁRIA. Sistema Brasileiro de Classificação de Solos. 2.ed. Rio de Janeiro: CNPSO, 2006. 306 p.

GOLDEN SOFTWARE. Surfer: user’s guide. Colorado: Golden Software, 2002. 664 p.

KANEVSKI, M.; MAIGNAN, M. Analysis and modelling of spatial environmental data. Lausanne: EPFL Press, 2004. 288 p.

LEAL, F. A.; MIGUEL, E. P.; MATRICARDI, E. A. T. Mapeamento de unidades produtivas utilizando a interpolação geoespacial krigagem a partir do inventário florestal em um povoamento de Eucalyptus urophylla S. T. Blake. Enciclopédia Biosfera, v. 7, n. 13, p. 727–745, 2011.

MATHERON, G. The theory of regionalized variables and its applications. Fontainebleau: École Nationale Supérieure des Mines de Paris, 1971. 211 p.

MELLO, J. M.; BATISTA, J. L. F.; OLIVEIRA, M. S.; RIBEIRO JR., P. J. Estudo da dependência espacial de características dendrométricas para Eucalyptus grandis. Cerne, Lavras, v. 11, n. 2, p. 113–126, 2005a.

MELLO, J. M.; BATISTA, J. L. F.; RIBEIRO JR., P. J.; OLIVEIRA, M. S. Ajuste e seleção de modelos espaciais de semivariograma visando à estimativa volumétrica de Eucalyptus grandis. Scientia Forestalis, Piracicaba, n. 69, p. 25–37, 2005b.

MOTOMIYA, A. V. A.; CORÁ, J. E.; PEREIRA, G. T. Uso da krigagem indicatriz na avaliação de indicadores de fertilidade do solo. Revista Brasileira de Ciência do Solo, Viçosa, v. 30, n. 3, p. 485–496, 2006.

NANOS, N.; CALAMA, R.; MONTERO, G.; GIL, L. Geostatistical prediction of height/diameter models. Forest Ecology and Management, Amsterdan, v. 195, n. 1-2, p. 221–235, 2004.

NANOS, N.; PARDO, F.; NAGER, J. A.; PARDOS, A. J.; GIL, L. Using multivariate factorial kriging for multiscale ordination: a case study. Canadian Journal of Forest Research, Ottawa, v. 35, n. 12, p. 2860-2874, 2005.

PAZ-FERREIRO, J.; VÁZQUEZ, E. V.; VIEIRA, S. R. Geostatistical analysis of a geochemical dataset. Bragantia, Campinas, v. 69, p. 121–129, 2010.

PELISSARI, A. L.; CALDEIRA, S. F.; DRESCHER, R.; SANTOS, V. S. Modelagem geoestatística da dinâmica espacial da altura dominante de Tectona grandis L.f. (teca). Enciclopédia Biosfera, v. 8, n. 15, p. 1249–1260, 2012.

PEREIRA, J. C.; MOURÃO, D. A. C.; SCALET, V.; SOUZA, C. A. M. Comparação entre modelos de relação hipsométrica com e sem componente espacial para Pinus sp. na FLONA Ipanema – SP. Scientia Forestalis, Piracicaba, v. 39, n. 89, p. 43–52, 2011.

REICHERT, J. M.; DARIVA, T. A.; REINERT, D. J.; SILVA, V. R. Variabilidade espacial de Planossolo e produtividade de soja em várzea sistematizada: análise geoestatística e análise de regressão. Ciência Rural, Santa Maria, v. 38, n. 4, p. 981–988, 2008.

RICHARDS, F. J. A. Flexible growth function for empirical use. Journal of Experimental Botany, Oxford, v. 10, n. 2, p. 290–300, 1959.

ROSA FILHO, G.; CARVALHO, M. P.; MONTANARI, R.; SILVA, J. M.; SIQUEIRA, G. M.; ZAMBIANCO, E. C. Variabilidade espacial de propriedades dendrométricas do eucalipto e de atributos físicos de um Latossolo Vermelho. Bragantia, Campinas, v. 70, n. 2, p.439–446, 2011.

RUFINO, T. M. C.; THIERSCH, C. R.; FERREIRA, S. O.; KANEGAE JR., H.; FAIS, D. Uso da Geoestatística no estudo da relação entre variáveis dentrométricas de povoamentos de Eucalyptus sp. e atributos do solo. Ambiência, Guarapuava, v. 2, n. 1, p. 83–93, 2006.

SCOLFORO, J. R. S. Biometria florestal: modelos de crescimento e produção florestal. Lavras: UFLA/FAEPE, 2006. 393 p.

SPURR, S. H. Forest inventory. New York: The Ronald Press Company, 1952. 476 p.

VIEIRA, S. R. Uso de geoestatística em estudos de variabilidade espacial de propriedades do solo. In: NOVAIS, R. F.; ALVAREZ V., V. H.; SCHAEFER, C. E. G. R. (Org.). Tópicos em ciência do solo. Viçosa: Sociedade Brasileira de Ciência do Solo, 2000. v. 1, p. 1–54.

WEBSTER, R.; OLIVER, M. A. Geostatistics for environmental scientists. 2.ed. West Sussex: John Wiley & Sons Ltd, 2007. 333 p.

YAMAMOTO, J. K.; LANDIM, P. M. B. Geoestatística: conceitos e aplicações. São Paulo: Oficina de textos, 2013. 215 p.