Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-1
5
GEOESTATíSTICA: FUNDAMENTOS E APLICAÇÕES
Eduardo Celso Gerbi Camargo
5.1 INTRODUÇÃO
Como nos ensina Burrough (1998), quando os dados são abundantes, a maior parte
dos métodos de interpolação produz valores semelhantes. Como discutido no capítulo
anterior, os métodos tradicionais de interpolação espacial, como triangulação, média
local das amostras e método da distância inversa, estão amplamente disponíveis nos
SIGs do mercado. No caso de dados esparsos, no entanto, tais métodos possuem
limitações na representação da variabilidade espacial, porque desconsideram a
anisotropia e a continuidade do fenômeno que se quer observar. Além disso, deixam
sem resposta algumas questões importantes, tais como (Burrough, 1987):
· Qual o tamanho ideal do domínio ou da janela de estimação ?
· Que forma e orientação deve ter a janela para se obter uma estimação ótima ?
· Existem outros modos para estimar os pesos além daqueles baseados em função
de distância ?
· Quais são os erros (incertezas) associados aos valores estimados ?
Na realidade, as propriedades naturais da superfície terrestre são espacialmente
contínuas, sendo restritivo descrevê-las através de simples funções matemáticas que não
respondem as questões acima formuladas. Modelos inferenciais para este objetivo vêm
sendo propostos. A krigeagem é um desses modelos, e sua base conceitual está
fundamentada na teoria das variáveis regionalizadas, formalizada por Matheron (1971).
O termo krigeagem é derivado do nome de Daniel G. Krige, que foi o pioneiro
em introduzir o uso de médias móveis para evitar a superestimação sistemática de
reservas em mineração (Delfiner e Delhomme, 1975). O que diferencia a krigeagem de
outros métodos de interpolação é a estimação de uma matriz de covariância espacial que
determina os pesos atribuídos às diferentes amostras, o tratamento da redundância dos
dados, a vizinhança a ser considerada no procedimento inferencial e o erro associado ao
valor estimado. Além disso, a krigeagem também fornece estimadores exatos com
propriedades de não tendenciosidade e eficiência.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-2
Portanto, a incorporação de procedimentos geoestatísticos em SIG’s, baseados
em técnicas de krigeagem, é importante, porque essa associação melhora os
procedimentos tradicionais de tais sistemas devido à qualidade do estimador e,
principalmente, pela informação de acurácia fornecida nesse modelo inferencial.
Neste contexto, este capítulo discute o uso de procedimentos geoestatísticos
baseado em técnicas de krigeagem no ambiente de Geoprocessamento. Estes
procedimentos englobam análise exploratória, geração e modelagem de
semivariogramas, validação do modelo e interpolação de krigeagem e tem como
principal objetivo a análise em duas dimensões para dados espacialmente distribuídos,
no que diz respeito à interpolação de superfícies geradas a partir de amostras
georreferenciadas obtidas por amostragem.
Os resultados aqui discutidos serviram de base para o desenvolvimento de um
Módulo de Procedimentos Geoestatísticos no SPRING (Camargo, 1997). Ferramentas
semelhantes estão disponíveis no ARC/INFO e GRASS.
5.2 FUNDAMENTAÇÃO TEÓRICA
A variabilidade espacial de algumas características do solo vem sendo uma das
preocupações de pesquisadores praticamente desde o início do século. Smith (1910)
estudou a disposição de parcelas no campo em experimentos de rendimento de
variedades de milho, numa tentativa de eliminar o efeito de variações do solo.
Montgomery (1913), preocupado com o efeito do nitrogênio no rendimento do trigo, fez
um experimento em 224 parcelas, medindo o rendimento de grãos. Vários outros
autores, como Waynick e Sharp (1919), também estudaram variações de nitrogênio e o
carbono no solo.
Os procedimentos usados na época baseavam-se na estatística clássica e
utilizavam grandes quantidades de dados amostrais, visando caracterizar ou descrever a
distribuição espacial da característica em estudo. Por estatística clássica entende-se
aquela que se utiliza de parâmetros como média e desvio padrão para representar um
fenômeno e se baseia na hipótese principal de que as variações de um local para outro
são aleatórias.
Krige (1951), trabalhando com dados de concentração de ouro, concluiu que
somente a informação dada pela variância seria insuficiente para explicar o fenômeno
em estudo. Para tal, seria necessário levar em consideração a distância entre as
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-3
observações. A partir daí surge o conceito da geoestatística, que leva em consideração a
localização geográfica e a dependência espacial.
Matheron (1963, 1971), baseado nas observações de Krige, desenvolveu a teoria
das variáveis regionalizadas. Uma variável regionalizada é uma função numérica com
distribuição espacial, que varia de um ponto a outro com continuidade aparente, mas
cujas variações não podem ser representadas por uma função matemática simples.
A teoria das variáveis regionalizadas pressupõe que a variação de uma variável
pode ser expressa pela soma de três componentes (Burrough, 1987): a) uma componente
estrutural, associada a um valor médio constante ou a uma tendência constante; b) uma
componente aleatória, espacialmente correlacionada; e c) um ruído aleatório ou erro
residual.
Se x representa uma posição em uma, duas ou três dimensões, então o valor da
variável Z, em x, é dada por (Burrough, 1987):
Z(x) = m(x) + e¢ (x) + e² (5.1)
onde:
· m(x) é uma função determinística que descreve a componente estrutural de Z
em x;
· e¢(x) é um termo estocástico, que varia localmente e depende espacialmente de
m(x);
· e² é um ruído aleatório não correlacionado, com distribuição normal com
média zero e variância s2.
As Figuras 5.1(a) e 5.1(b) ilustram as três componentes principais da variação
espacial. A Figura 5.1(a) apresenta uma componente determinística que varia
abruptamente, enquanto a componente determinística na Figura 5.1(b) apresenta uma
tendência constante.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-4
Figs. 5.1(a) e 5.1(b) - Principais componentes da variação espacial.
FONTE: Modificada de Burrough (1987), p. 155.
5.2.1 HIPÓTESES CONSIDERADAS
Diferente dos métodos convencionais de estimação, a krigeagem está
fundamentada na teoria das variáveis regionalizadas. A hipótese mais comum é a
chamada "estacionaridade de 2a. ordem" (Burrough, 1987):
· A componente determinística, m(x), é constante (não há tendências na região).
· A variância das diferenças entre duas amostras depende somente da distância
h entre elas, isto é:
Var[Z(x)-Z(x+h)] = E{[Z(x)-Z(x+h)]2} = 2g(h) , (5.2)
onde g(h) é chamado de semi-variância.
Para mostrar a contribuição da semivariância, podemos reescrever a equação 5.1 como
Z(x) = m(x) + g(h) + e² (5.3)
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-5
Em outras palavras, como supomos m(x) ser constante, a variação local das amostras (e
seu relacionamento espacial) pode ser caracterizado pela semivariância g(h).
5.3 VARIOGRAMA
O variograma é uma ferramenta básica de suporte às técnicas de krigeagem, que permite
representar quantitativamente a variação de um fenômeno regionalizado no espaço
(Huijbregts, 1975). Considere duas variáveis regionalizadas, X e Y, onde X = Z(x) e Y
= Z(x+h). Neste caso, referem-se ao mesmo atributo (por exemplo, o teor de zinco no
solo) medido em duas posições diferentes, conforme ilustra a Figura 5.2, onde
Fig. 5.2 - Amostragem em duas dimensões.
x denota uma posição em duas dimensões, com componentes (xi , yi), e h um vetor
distância (módulo e direção) que separa os pontos.
O nível de dependência entre essas duas variáveis regionalizadas, X e Y, é
representado pelo variograma, 2g(h), o qual é definido como a esperança matemática do
quadrado da diferença entre os valores de pontos no espaço, separados pelo vetor
distância h, isto é,
2g(h) = E{[Z(x)-Z(x+h)]2} = Var[Z(x)-Z(x+h)] . (5.4)
Através de uma amostra z(xi), i=1, 2, ..., n, o variograma pode ser estimado por
2 ( =
1
N( )
z( )- z( )
N( )
g 2
Ù
h)
h
x x h
h
[ i i ]
i=
å +
1
, (5.5)
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-6
onde:
· 2g
Ù
(h) - é o variograma estimado;
· N(h) - é o número de pares de valores medidos, z(xi) e z(xi+h), separados por
um vetor distância h;
· z(xi) e z(xi+h), - são valores da i-ésima observação da variável
regionalizada, coletados nos pontos xi e xi+h (i = 1, ..., n), separados pelo vetor
h.
5.3.1 PARÂMETROS DO SEMIVARIOGRAMA
A Figura 5.3 ilustra um semivariograma experimental com características muito
próximas do ideal. O seu padrão representa o que, intuitivamente, se espera de dados de
campo, isto é, que as diferenças {Z(xi) - Z(xi + h)} decresçam à medida que h, a
distância que os separa decresce. É esperado que observações mais próximas
geograficamente tenham um comportamento mais semelhante entre si do que aquelas
separadas por maiores distâncias. Desta maneira, é esperado que g(h) aumente com a
distância h.
Fig. 5.3 - Exemplo de semivariograma.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-7
Os parâmetros do semivariograma podem ser observados diretamente da Figura
5.3 :
· Alcance (a): distância dentro da qual as amostras apresentam-se
correlacionadas espacialmente. Na Figura 5.3, o alcance ocorre próximo de
25m.
· Patamar (C): é o valor do semivariograma correspondente a seu alcance (a).
Deste ponto em diante, considera-se que não existe mais dependência espacial
entre as amostras, porque a variância da diferença entre pares de amostras
(Var[Z(x) - Z(x+h)]) torna-se invariante com a distância.
· Efeito Pepita (C0): idealmente, g(0)=0. Entretanto, na prática, à medida que h
tende para 0 (zero), g(h) se aproxima de um valor positivo chamado Efeito
Pepita (C0), que revela a descontinuidade do semivariograma para distâncias
menores do que a menor distância entre as amostras. Parte desta
descontinuidade pode ser também devida a erros de medição (Isaaks e
Srivastava, 1989), mas é impossível quantificar se a maior contribuição provém
dos erros de medição ou da variabilidade de pequena escala não captada pela
amostragem.
· Contribuição (C1): é a diferença entre o patamar (C) e o Efeito Pepita (Co).
5.3.2 CÁLCULO DO SEMIVARIOGRAMA A PARTIR DE AMOSTRAS REGULARMENTE
ESPAÇADAS
Considere o conjunto de amostras regularmente espaçadas, em duas dimensões,
conforme apresentado na Figura 5.4.
Fig. 5.4 - Amostras regularmente espaçadas em duas dimensões.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-8
Para determinar o semivariograma experimental, por exemplo, na direção de 900
o cálculo de g
Ù
(h)é repetido para todos os intervalos de h. Suponha a distância entre
dois pontos consecutivos igual a 100 metros (d=100m). Então, qualquer par de
observações, na direção 900, cuja distância é igual a 100m será incluído no cálculo de
gÙ
(90o , 100m) . Isto feito, os cálculos são repetidos para a próxima distância, por
exemplo, 200m. Isto inclui todos os pares de observações cuja distância é igual a 200m.
O processo é repetido até que algum ponto de parada desejado seja alcançado. Este
procedimento pode ser melhor compreendido com o auxílio da Figura 5.5 e também
deve ser realizado para outras direções (00, 450 e 1350).
Fig.5.5 - Ilustração para o cálculo do semivariograma a partir de amostras regularmente
espaçadas.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-9
5.3.3 CÁLCULO DO SEMIVARIOGRAMA A PARTIR DE AMOSTRAS
IRREGULARMENTE ESPAÇADAS
Considere o conjunto de amostras irregularmente espaçadas, em duas
dimensões, conforme apresentado na Figura 5.6. Neste caso, para determinar o
semivariograma experimental, é necessário introduzir limites de tolerância para direção
e distância.
Fig. 5.6 - Parâmetros para o cálculo do semivariograma a partir de amostras
irregularmente espaçadas em duas dimensões.
FONTE: Modificada de Deutsch e Journel (1992)
Tome como referência o Lag2 (Lag refere-se a uma distância pré-definida, a
qual é utilizada no cáculo do semivariograma) da figura acima. Suponha um
incremento de Lag igual a 100 metros com tolerância de 50 metros. Considere ainda a
direção de medida 450 com tolerância angular 22.50. Então, qualquer par de observações
cuja distância está compreendida entre 150m e 250m e 22.50 e 67.50 será incluído no
cálculo do semivariograma de Lag2. Este processo se repete para todos os Lag’s.
Ainda com referência na Figura 5.6, a largura de banda (BW) se refere a um
valor de ajuste a partir do qual se restringe o número de pares de observações para o
cálculo do semivariograma. A próxima etapa constitui o ajuste de um modelo teórico ao
semivariograma experimental, conforme descrito a seguir.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-10
5.3.4 MODELOS TEÓRICOS
O gráfico do semivariograma experimental, g
Ù
(h), calculado através da Equação
(5.5), é formado por uma série de valores, conforme ilustra a Figura 5.3, sobre os quais
se objetiva ajustar uma função. É importante que o modelo ajustado represente a
tendência de g
Ù
(h) em relação a h. Deste modo, as estimativas obtidas a partir da
krigeagem serão mais exatas e, portanto mais confiáveis.
O procedimento de ajuste não é direto e automático, como no caso de uma
regressão, por exemplo, mas sim interativo, pois nesse processo o intérprete faz um
primeiro ajuste e verifica a adequação do modelo teórico. Dependendo do ajuste obtido,
pode ou não redefinir o modelo, até obter um que seja considerado satisfatório.
Os modelos aqui apresentados são considerados modelos básicos, denominados
de modelos isotrópicos por Isaaks e Srivastava (1989). Estão divididos em dois tipos:
modelos com patamar e modelos sem patamar. Modelos do primeiro tipo são
referenciados na geoestatística como modelos transitivos. Alguns dos modelos
transitivos atingem o patamar (C) assintoticamente. Para tais modelos, o alcance (a) é
arbitrariamente definido como a distância correspondente a 95% do patamar. Modelos
do segundo tipo não atingem o patamar, e continuam aumentanto enquanto a distância
aumenta. Tais modelos são utilizados para modelar fenômenos que possuem
capacidade infinita de dispersão. Os modelos transitivos mais utilizados são: modelo
esférico (Sph), modelo exponencial (Exp) e modelo gaussiano (Gau). Estes modelos
estão apresentados na Figura 5.7 com o mesmo alcance (a).
Fig. 5.7 - Representação gráfica de modelos transitivos normalizados.
FONTE: Modificada de Isaaks e Srivastava (1989).
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-11
5.3.4.1 MODELO ESFÉRICO
O modelo esférico é um dos modelos mais utilizados e está representado em vermelho
na Figura 5.7. A equação normalizada deste modelo é:
Sph( h )
h
h h
h|
h|
=
=
æ
è ç
ö
ø ÷
-
æ
è ç
ö
ø ÷
< £
>
ì
í
ïï
î
ïï
0
5 5 0
1
,
, , , |
, |
| | 0
1
| |
0
| | 3
a a
a
a
(5.6)
5.3.4.2 MODELO EXPONENCIAL
Um outro modelo bastante utilizado é o modelo exponencial, o qual é
apresentado em azul na Figura 5.7. A equação normalizada deste modelo é:
Exp( )
0 , |
1 exp
h | |
h
h
h
=
- -
æ
è ç
ö
ø ÷
¹
ì
í ï
î ï
|=0
|
a
, | 0
(5.7)
Este modelo atinge o patamar assintoticamente, com o alcance prático definido
como a distância na qual o valor do modelo é 95% do patamar (Isaaks e Srivastava,
1989).
5.3.4.3 MODELO GAUSSIANO
O modelo gaussiano é um modelo transitivo, muitas vezes usado para modelar
fenômenos extremamente contínuos (Isaaks e Srivastava, 1989). Sua formulação é dada
por:
Gau ( )
|=0
1 exp
| | 2
| 0
h
h
h
h
=
- -
æ
è ç
ö
ø ÷
¹
ì
í ï
î ï
0 , |
, |
a
(5.8)
Semelhante no modelo exponencial, o modelo gaussiano atinge o patamar
assintoticamente e o parâmetro a é definido como o alcance prático ou distância na qual
o valor do modelo é 95% do patamar (Isaaks e Srivastava, 1989). O que caracteriza este
modelo é seu comportamento parabólico próximo à origem, conforme representado na
Figura 5.7 através da linha sólida verde.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-12
5.3.4.4 MODELO POTÊNCIA
O modelo potência não é um modelo transitivo, portanto não atinge o patamar.
Em geral, este tipo de modelo é utilizado para modelar fenômenos com capacidade
infinita de dispersão. A Figura 5.8 ilustra o modelo potência, o qual é expresso através
de:
Pot (|h|)
h
h h
=
¹
ìí ï
î ï
0 ,|
| ,|
|=0
c. |e | 0 (5.9)
onde,
· c é o coeficiente de declividade, e
· e é o expoente.
Fig.5.8 - Representação gráfica do modelo potência.
Até este ponto foram apresentados os principais modelos básicos normalizados,
os quais são utilizados para modelar ou ajustar o semivariograma experimental. Na
prática, os semivariogramas experimentais possuem valores de efeito pepita (Co) maior
que zero e valores de patamar (C) maiores que a unidade, conforme ilustrado na Figura
5.9.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-13
Fig. 5.9 - Representação gráfica de semivariogramas experimentais e modelos teóricos.
Em resumo, os semivariogramas dos modelos transitivos básicos são assim definidos:
· Modelo Esférico de Semivariograma:
g( )
0 , | |= 0
3
2
| | 1
2
| | 3
Sph (| |) |
|
h
h
h h
h h
h
= +
æ
è ç
ö
ø ÷
-
æ
è ç
ö
ø ÷
é
ë êê
ù
û úú
= + < £
+ >
ì
í
ïï
î
ïï
Co C1 a a
Co C1 a
Co C1 a
[ ] , |
, |
0 (5.10)
· Modelo Exponencial de Semivariograma:
g( )
| | 0
1 exp
| |
Exp (| |)] ,| | 0
h
h
h
h h
=
=
- -
æ
è ç
ö
ø ÷
é
ë ê
ù
û ú
= ¹
ì
í ï
î ï
0 ,
Co+C1 [
a Co+C1
(5.11)
· Modelo Gaussiano de Semivariograma:
g( )
| | 0
1 exp
| | 2
[Gau (| |)] ,| | 0
h
h
h
h h
=
=
- -
æ
è ç
ö
ø ÷
é
ë êê
ù
û úú
= ¹
ì
í ï
î ï
0 ,
Co+C1
a Co+C1
(5.12)
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-14
De maneira análoga, o modelo potência é escrito em termos de semivariograma da
seguinte forma:
· Modelo Potência de Semivariograma :
g(h)
h
h h h
=
= ¹
ìí ï
î ï
0 ,|
| ,|
|=0
C0+c. |e C0+ Pot (| |) | 0
(5.13)
5.3.4.5 MODELOS ANINHADOS
Existem determinados fenômenos em que são necessários modelos mais
complexos de semivariograma para explicar suas variações espaciais. Estes modelos são
combinações de modelos simples, denominados aninhados. McBratney e Webster
(1986) observaram que modelos aninhados são necessários para explicar a variação do
solo decorrente de fatores independentes de formação. Por exemplo, um modelo
aninhado útil em estudos de mineração e pesquisa de solo é o duplo esférico.
McBratney et al. (1982) o utilizaram para descrever a variação do cobre e do cobalto no
solo. Este modelo é definido como:
g
g
g
1
2
( )
3
2
| | 1
2
3
| |
( ) | |
3
2
| | 1
2
3
| |
( ) | |
| |
| |=0
h
h h
h h
h h
h h
h
h
=
+
æ
è
çç
ö
ø
÷÷
-
æ
è
çç
ö
ø
÷÷
é
ë
êê
ù
û
úú
= < £
+
æ
è
çç
ö
ø
÷÷
-
æ
è
çç
ö
ø
÷÷
é
ë
êê
ù
û
úú
= < £
+ >
ì
í
ïïïïï
î
ïïïïï
C C
a a
a
C C
a a
a a
C C a
0 1
1 1
1
0 2
2 2
1 2
0 2 2
,
,
,
,
0
0
(5.14)
onde,
· a1 e C1 correspondem aos parâmetros de alcance e contribuição,
respectivamente, do primeiro modelo esférico (g1 (h)) e
· a2 e C2 correspondem aos parâmetros de alcance e contribuição,
respectivamente, do segundo modelo esférico (g2 (h)).
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-15
Este modelo é mostrado na Figura 5.10, onde as linhas sólidas representam os
modelos de ajuste teórico ao semivariograma experimental.
Fig. 5.10 - Representação gráfica de um modelo duplo esférico.
Dependendo do fenômeno em estudo, outros modelos aninhados são necessários
para caracterizar a variabilidade espacial. Por exemplo: duplo exponencial, exponencial
com duplo esférico, linear com duplo esférico, etc.
5.3.5 ANISOTROPIA
A anisotropia pode ser facilmente constatada através da observação dos
semivariogramas obtidos para diferentes direções. As convenções direcionais usadas na
geoestatística são mostradas na Figura 5.11.
Fig. 5.11 - Convenções direcionais usadas na geoestatística.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-16
Considere os semivariogramas obtidos para as direções 00, 450, 900 e 1350,
ilustrados na Figura 5.15. Verifica-se uma similaridade bastante grande entre eles. Esta
é a representação de um caso simples e menos freqüente, em que a distribuição espacial
do fenômeno é denominada isotrópica. Neste caso, um único modelo é suficiente para
descrever a variabilidade espacial do fenômeno em estudo.
Fig. 5.12 - Representação gráfica de semivariogramas isotrópicos.
Por outro lado, se os semivariogramas não são iguais em todas as direções, a
distribuição é denominada anisotrópica. Se a anisotropia é observada e é refletida pelo
mesmo Patamar (C) com diferentes Alcances (a) do mesmo modelo, então ela é
denominada Geométrica.
Considere o semivariograma ilustrado na Figura 5.13. Os pontos interligados
com linhas tracejadas são os semivariogramas experimentais em duas direções
ortogonais. O semivariograma que atinge primeiro o patamar (azul) se refere à direção
de 1200 e o semivariograma com maior alcance (vermelho) se refere à direção de 300.
As linhas sólidas em ambas direções são os modelos teóricos de ajuste dos
semivariogramas experimentais.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-17
Fig. 5.13 - Representação gráfica de anisotropia geométrica.
Um modo direto de visualizar e calcular os parâmetros (fator e ângulo) da
anisotropia geométrica é através do esboço gráfico de uma elipse, calculada através dos
alcances obtidos em direções distintas, conforme Figura 5.14. As convenções que
seguem, são as adotadas por Deutsch e Journel (1992). Para o eixo maior da elipse,
denominado direção de máxima continuidade, aplica-se o maior alcance(a1). O ângulo
da direção de máxima continuidade é definido a partir da direção Norte e no sentido
horário. Seu valor corresponde à direção de maior alcance. O eixo menor define o
alcance(a2) na direção de menor continuidade, sendo este ortogonal à direção principal.
Fig. 5.14 - Representação gráfica da anisotropia geométrica em duas dimensões.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-18
FONTE: Modificada de Deutsch e Journel (1992), p. 24.
O fator de anisotropia geométrica é definido como a razão entre o alcance na
direção de menor continuidade (a2) e o alcance na direção de maior continuidade (a1).
Neste caso, o fator de anisotropia geométrica é sempre menor que a unidade e o ângulo
de anisotropia é igual ao ângulo da direção de máxima continuidade.
Existe ainda um outro tipo de anisotropia em que os semivariogramas
apresentam os mesmos Alcances (a) e diferentes Patamares (C). Neste caso, a
anisotropia é denominada Zonal. Como a isotropia, a anisotropia zonal também é um
caso menos freqüente presente nos fenômenos naturais. O mais comum é encontrar
combinações da anisotropia zonal e geométrica, denominada anisotropia combinada.
Considere o semivariograma apresentado na Figura 5.15. Os pontos interligados
com linhas tracejadas correspondem a semivariogramas experimentais em duas direções
ortogonais. O semivariograma com maior patamar (azul) refere-se à direção de 600 e o
semivariograma com menor patamar (vermelho) refere-se à sua direção perpendicular
(1500). Os modelos de ajuste aos semivariogramas estão representados por linhas
sólidas.
Fig. 5.15 - Representação gráfica de anisotropia combinada.
Segundo Isaaks e Srivastava (1989), a anisotropia zonal pode ser considerada
como um caso particular da anisotropia geométrica, ao se supor um fator de anisotropia
muito grande. Nesta condição, o alcance implícito na direção de menor continuidade é
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-19
muito grande. A estrutura do semivariograma é então adicionada somente para a direção
de maior continuidade.
5.4 KRIGEAGEM
O termo krigeagem é derivado do nome Daniel G. Krige, que foi o pioneiro a
introduzir o uso de médias móveis para evitar a superestimação sistemática de reservas
de mineração (Delfiner e Delhomme, 1975). Inicialmente, o método de krigeagem foi
desenvolvido para solucionar problemas de mapeamentos geologicos, mas seu uso
expandiu-se com sucesso no mapeamento de solos (Burgess e Webster, 1980a,b),
mapeamento hidrológico (Kitanidis e Vomvoris, 1983), mapeamento atmosférico
(Lajaunie, 1984) e outros campos correlatos.
A diferença entre a krigeagem e outros métodos de interpolação é a maneira
como os pesos são atribuídos às diferentes amostras. No caso de interpolação linear
simples, por exemplo, os pesos são todos iguais a 1/N (N = número de amostras); na
interpolação baseada no inverso do quadrado das distâncias, os pesos são definidos
como o inverso do quadrado da distância que separa o valor interpolado dos valores
observados. Na krigeagem, o procedimento é semelhante ao de interpolação por média
móvel ponderada, exceto que aqui os pesos são determinados a partir de uma análise
espacial, baseada no semivariograma experimental. Além disso, a krigeagem fornece,
em média, estimativas não tendenciosas e com variância mínima1.
Segundo Oliver e Webster (1990), a krigeagem engloba um conjunto de métodos
de estimação, a saber: krigeagem simples, krigeagem ordinária, krigeagem universal,
Co-krigeagem, krigeagem disjuntiva, etc. Este trabalho limita-se somente à
apresentação da krigeagem ordinária, descrita a seguir.
1Estimativas não tendenciosas significam que, em média, a diferença entre valores
estimados e verdadeiros para o mesmo ponto deve ser nula; e variância mínima
significa que estes estimadores possuem a menor variância dentre todos os estimadores
não tendenciosos.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-20
5.4.1 KRIGEAGEM ORDINÁRIA
Considere uma superfície sobre a qual se observe alguma propriedade do solo, Z,
em n pontos distintos, com coordenadas representadas pelo vetor x. Assim, tem-se um
conjunto de valores {z(xi), i=1, ..., n}, onde xi, identifica uma posição em duas
dimensões representada pelos pares de coordenadas (xi, yi). Suponha que se objetive
estimar o valor de Z no ponto x0. O valor desconhecido de Z(x0) pode ser estimado a
partir de uma combinação linear dos n valores observados, adicionado a um parâmetro,
l0 (Journel, 1988):
Z* Z( )
0 0 i x x = +å
=
l li
n
i 1
. (5.15)
Deseja-se um estimador não tendencioso, isto é,
E [Z xo- Z xo
* ] = 0 . (5.16)
A relação acima impõe que as duas médias sejam iguais; assim aplicando-se a
Equação 5.39 em 5.40, obtem-se:
E [Z ] E i
i
n
i
i
n
x x o = + Z( i
é
ë ê
ù
û ú
Þ = +
= =
l å l l å l 0
1
0
1
. ) m m . (5.17)
A krigeagem ordinária não requer o prévio conhecimento da média m. Neste caso,
para que a igualdade da Equação 5.41 seja satisfeita é necessário que:
l l 0 = 0 e i=
= å
i
n
1
1 .
Portanto, o estimador de krigeagem ordinária é:
Z*
0 x x = å
=
li
n
i
i
Z
1
( ) , com li
i
n
= å
=
1
1 . (5.18)
Journel (1988) mostra que, minimizando a variância do erro (Var[Z Z* ]
x x 0 0
- ) sob a
condição de que li
i
n
= å
=
1
1 , os pesos li são obtidos a partir do seguinte sistema de
equações, denominado sistema de krigeagem ordinária:
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-21
l a
l
j i j i 0
j
C C( para i 1,..., n
j 1
n
j 1
n
( , ) , )
=
=
å
å
- = =
=
ì
í
ïïï
î
ïïï
x x x x
1
(5.19)
onde,
· C(xi, xj) e C(xi, x0) são respectivamente, a semivariância entre os pontos xi e xj
e entre os pontos xi e x0.
· a é o multiplicador de Lagrange necessário para a minimização da variância
do erro.
A correspondente variância minimizada do erro, denominada variância de krigeagem
ordinária (sko
2 ), é dada pela seguinte expressão (Journel, 1988):
s l a ko
2
i
i=1
n
= Var [Z - Z* ]= C( ) -å C( , ) - x x i 0 x x
0 0 0. (5.20)
A krigeagem ordinária é um interpolador exato no sentido de que, quando as equações
acima forem usadas, os valores interpolados irão coincidir com os valores dos pontos
amostrais. A variância da krigeagem ordinária, indicada na equação (5.20), fornece
informação importante sobre a confiabilidade dos valores interpolados (Burrough,
1998).
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-22
5.5 EXEMPLO DE APLICAÇÃO
Esta seção mostra, através de um experimento prático, a importância da integração de
um SIG com procedimentos geoestatísticos na análise e interpolação de dados
espacialmente distribuídos.
Os dados utilizados, de propriedade do Centro Nacional de Pesquisas de Solos (CNPS -
RJ), foram obtidos no levantamento dos solos da Fazenda Canchim, em São Carlos -
SP. Estes se referem a uma amostragem de 85 observações georreferenciadas coletadas
no horizonte Bw (camada do solo com profundidade média de 1m). Dentre as variáveis
disponíveis, selecionou-se para estudo o teor de argila.
O estudo de caso aqui apresentado refere-se à análise da variação espacial do teor de
argila sobre a área da Fazenda Canchim. Considera-se o teor de argila ao longo do
perfil, classificado do seguinte modo (Calderano Filho et al., 1996):
· MUITO ARGILOSO: solos que apresentam 59% ou mais de argila;
· ARGILOSO: solos que apresentam de 35% a 59% de argila;
· MÉDIO: solos que apresentam de 15% a 35% de argila;
· ARENOSO: solos que apresentam menos de 15% de argila.
Para tal experimento, foram realizados os procedimentos de análise geoestatística:
análise exploratória e de continuidade espacial do teor de argila, modelagem dos
semivariogramas experimentais, validação do modelo e interpolação através de
krigeagem. Para maiores detalhes, o leitor deve referir-se a Camargo (1997), que
descreve ainda o desenvolvimento de um módulo geoestatístico no ambiente SPRING.
5.5.1 ANÁLISE EXPLORATÓRIA
Como mencionado anteriormente, a amostra é composta por 85 observações
georreferenciadas; dentre as variáveis disponíveis, escolheu-se o teor de argila. O valor
desta propriedade corresponde às observações que se realizaram no interior do solo,
horizonte Bw. A Figura 5.16 apresenta a distribuição das observações na área de
estudo.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-23
Fig 5.16 - Distribuição das observações na área de estudo.
Utilizando os recursos disponíveis do Módulo de Análise Exploratória desenvolvido no
sistema SPRING, pode-se determinar que a distribuição da variável em estudo é
aproximadamente normal, dando suporte à hipótese de normalidade dos dados.
5.5.2 MODELAGEM DOS SEMIVARIOGRAMAS E DA ANISOTROPIA
Após análise exaustiva das diferente configurações dos dados, descrita em Barbosa
(1997), foram considerados dois semivariogramas, como representativos da anisotropia
dos dodos: o associados à direção Norte (0o), que possui maior alcance e o
semivariograma associado à direção Leste (90o), que possui menor alcance. O resultado
do ajuste experimental dos semivariogramas é apresentado na Figura 5.17.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-24
Fig. 5.17 - Representação gráfica dos modelos teóricos de ajuste aos semivariogramas
experimentais nas direções 0o e 90o.
Conforme ilustrado na Figura 5.17, as linhas azuis, em ambas direções, são os modelos
teóricos de ajuste dos semivariogramas direcionais. O semivariograma com maior
patamar foi ajustado com um efeito pepita de (83,0), uma estrutura esférica com
contribuição (352,0) e alcance (4.031,0). De maneira análoga, o semivariograma na
direção 90o foi ajustado com efeito pepita (45,5), uma estrutura esférica com
contribuição (210,0) e alcance (1.883,0).
O passo seguinte refere-se a modelagem da anisotropia. Inicialmente identifica-se os
eixos de anisotropia, isto é, os eixos de maior e menor variabilidade espacial do
fenômeno em estudo. Neste caso, em particular, os eixos de anisotropia coincidem com
os eixos de coordenadas. Uma vez identificados os eixos de anisotropia, define-se um
único modelo que descreve como o semivariograma varia quando mudam as distâncias
e direções.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-25
Um método para combinar os modelos direcionais num único modelo consistente para
todas as direções é apresentado por Isaaks e Srivastava (1989, p. 377-386).
Resumidamente, este método consiste em definir uma transformação que reduz o
modelo direcional a um modelo normalizado com alcance igual a 1 (substitui-se h pela
distância h/a, onde a é o alcance). O modelo único (ou completo), o qual é expresso
através da seguinte equação:
g 4031,0 4031,0
( ) , , ,
,0
h = + + + + +
æ
è ç
ö
ø ÷
æ
è
çç
ö
ø
÷÷
é
ë
êê
ù
û
úú
¥
æ
è
çç
ö
ø
÷÷
é
ë
êê
ù
û
úú
æ
è çç
ö
ø ÷÷
æ
è
çç
ö
ø
÷÷
æ
è
çç
ö
ø
÷÷
45 5 210 0
2 2
142 0
2 2
1883
Sph
h0o h90o
Sph
h0o h90o
+
æ
è ç
ö
ø ÷
+
¥
æ
è ç
ö
ø ÷
æ
è
çç
ö
ø
÷÷
é
ë
êê
ù
û
úú
37 5
2 2
, Sph h0o h90o
e
. (5.21)
O modelo completo, apresentado na Equação 4.3, é composto de quatro termos, a saber:
· O primeiro termo referente ao efeito pepita (Co), igual a 45,5.
· O segundo termo referente à 1a componente estrutural, um modelo esférico
com contribuição 210, alcance 4031,0 na direção Norte (0o) e alcance 1883,0
na direção Leste (90o).
· O terceiro termo referente à 2a componente estrutural, um modelo esférico
com contribuição 142,0 e alcance 4031,0 relativo à direção Norte (0o).
· Por último, o quarto referente à 3a componente estrutural, um modelo esférico
com contribuição 37,5 e alcance muito pequeno, e, relativo à direção Norte
(0o). Este é um artifício empregado para representar o segundo efeito pepita .
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-26
5.5.3 ANÁLISE E APRESENTAÇÃO DE RESULTADOS
Para a apresentação dos resulatdos, é adotada a nomenclatura abaixo descrita:
· GKrV: refere-se à grade de valores estimados obtida da interpolação de
krigeagem, a partir de um modelo anisotrópico que supostamente representa a
verdadeira continuidade espacial do fenômeno em estudo. Da análise
geoestatística realizada, supõe-se que o fenômeno apresenta maior
continuidade na direção Norte (0o) e menor na direção Leste (90o).
· GKrVe: refere-se à grade da variância de krigeagem, associada a GKrV.
· GKrI: refere-se a Grade de valores estimados obtida da interpolação de
krigeagem a partir de um modelo isotrópico. Neste caso, admite-se que a
continuidade espacial do fenômeno é a mesma em qualquer direção.
· GKrIe: refere-se à grade da variância de krigeagem, associada a GKrI.
· GKrA : refere-se à grade de valores estimados obtida da interpolação de
krigeagem, a partir de um modelo anisotrópico que utiliza direções
intermediárias (10o e 100o) às direções de máxima e mínima continuidade.
· GKrAe: refere-se à grade da variância de krigeagem, associada a GKrA.
5.5.4 RESULTADOS DA KRIGEAGEM
Como ponto de partida, é interessante verificar a continuidade espacial do
fenômeno em estudo (teor de argila). Para realizar tal análise, é necessário transformar
as grades de valores estimados e as correspondentes grades de erros em imagens.
A seqüência de imagens a seguir apresenta os resultados obtidos sobre as grades
GKrV (Figura 5.18), GKrVe (Figura 5.19), GKrI (Figura 5.20), GKrIe (Figura 5.21),
GKrA (Figura 5.22) e GKrAe (Figura 5.23).
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-27
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-28
Através das imagens apresentadas, pode-se constatar algumas características
comuns:
· essas imagens revelam que, nas regiões Norte e Nordeste, o teor de argila é
relativamente baixo;
· na região Central, observa-se mudanças graduais do teor de argila, indo de
valores moderados a altos, e
· nas regiões Sul e Sudoeste aproximadamente moderado.
De maneira análoga, as grades de variância do erro mostram que o erro da estimativa
aumenta à medida em que se afasta dos pontos de observações. É possível também
identificar, nestas imagens, regiões onde a amostragem pode ser melhorada.
Um outro aspecto a ser observado, e talvez o mais importante, é o efeito da
anisotropia:
· Observa-se, na grade GKrV, a qual supostamente representa a verdadeira
continuidade espacial, que as mudanças graduais do teor de argila são
visivelmente diferentes das demais, principalmente na região central.
· A imagem da grade GkrI mostra que a continuidade espacial do teor de argila se
propaga uniformemente em todas as direções. Neste caso a anisotropia é
mascarada e, portanto, o resultado não revela a verdadeira continuidade espacial
da variável em estudo.
· Por outro lado, a imagem da grade GkrA apresenta um caso intermediário à da
suposta continuidade espacial verdadeira. Este caso também não revela a
verdadeira variabilidade espacial, apenas mostra que o teor de argila se propaga
mais intensamente na direção 10o e menos intensamente na direção ortogonal
(100o).
Resumindo, esta primeira análise procurou enfocar, através de três exemplos, que a
escolha do modelo decorrente da análise de continuidade espacial deve ser a melhor
possível, a fim de retratar a verdadeira variabilidade do fenômeno em estudo.
A análise seguinte refere-se a uma operação de fatiamento realizada na grade GKrV, a
qual supostamente representa a verdadeira continuidade espacial do teor de argila. Os
valores utilizados no fatiamento estão definidos segundo Calderano Filho et al. (1996).
O resultado desta operação é a imagem apresentada na Figura 5.25.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-29
Fig. 5.24 - Geologia da região de estudo.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-30
Fig. 5.25 - Imagem relativa ao fatiamento
da grade GKrV.Observando as Figuras 5.24 e 5.25, nota-se que ambas
possuem feições semelhantes. Segundo Calderano Filho et al. (1996), no Arenito
Superficial, encontram-se os solos argilosos e de teor médio de argila, na Formação
Serra Geral (diabásio) predominam os solos com teores argilosos a muito argilosos e na
Formação Arenito Botucatu encontram-se solos com teor médio de argila e os solos
arenosos.
Um outro tipo de análise refere-se à comparação entre o mapa de classes de
solos segundo o teor de argila, Mp1, (concebido pelo CNPS-EMBRAPA) e o mapa de
teor de argila, Mp2, (gerado através de procedimentos geoestatísticos). Neste caso, a
análise é restrita, pois, considera-se apenas o teor de argila às respectivas classes de
solos. A Figura 5.26 confronta os resultados expressos pelos mapas Mp1 e Mp2.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-31
Fig. 5.26 - Comparação entre o mapa de classes de solos segundo o teor de argila e o
mapa de teor de argila gerado através de procedimentos geoestatísticos.
Na Figura 5.26, alguns resultados são observados. Por exemplo, examinando a
região norte e nordeste do mapa de solos (Mp1), observa-se a presença de solos
arenosos (classe AQd) e solos com teor médio de argila (classes LEd3, LEa2, LEa3), o
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-32
que pode ser constatado também no mapa de teor de argila (Mp2). Um outro exemplo,
são os solos argilosos e muito argilosos (classes de solos Lud, LRd1 e Tre1), que se
encontram na região central de ambos os mapas, Mp1 e Mp2. E assim, várias outras
classes de solo, como LVa1, LVa2, etc.
Considerando a deficiência de amostragem, com as conseqüentementes
limitações impostas à análise geoestatística efetuada, o resultado apresentado na Figura
5.26 é de certa forma expressivo e revela a boa adequação das técnicas de krigeagem ao
estudo de caso.
5.5.5 KRIGEAGEM VERSUS MÉTODOS CONVENCIONAIS DE ESTIMAÇÃO
Esta seção apresenta, de forma comparativa, os resultados obtidos da krigeagem
com outros métodos de estimação disponíveis no sistema SPRING, como: média
ponderada pelo inverso do quadrado da distância, média simples e vizinho mais
próximo. Não é o objetivo deste trabalho detalhar cada um destes métodos, pois os
mesmos são descritos na literatura, como, por exemplo, em Burrough (1987).
Das 85 observações de campo realizadas nos limites da Fazenda Canchim, as seguintes
grades de valores estimados são geradas:
· GKrV e GKrI : conforme definido anteriormente;
· GID: refere-se à grade de valores estimados obtida através do método de
estimação média ponderada pelo inverso do quadrado da distância;
· GMS: refere-se à grade de valores estimados obtida através do método de
estimação média simples;
· GVP: refere-se à grade de valores estimados obtida através do método de
estimação vizinho mais próximo.
A seqüência de imagens a seguir apresenta os resultados obtidos e se refere às grades
GKrV , GKrI , GID, GMS e GVP . O mapa geológico é utilizado como referência à
análise comparativa.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-33
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-34
Através das imagens apresentadas , pode-se constatar alguns
fatos:
· É relativamente fácil de ser constatado que a imagem referente à GKrV
(Grade de Krigeagem obtida a partir de um modelo anisotrópico e que
supostamente representa a Verdadeira continuidade espacial do fenômeno) é a
que melhor expressa a variabilidade do teor de argila. Por exemplo, no mapa
geológico onde ocorre o diabásio, cujo teor de argila é do produto de alteração
grande, observa-se que o resultado expresso pela grade GKrV tende a moldar
a variabilidade espacial do fenômeno nesta região.
· Com algumas ressalvas, o método da média ponderada pelo inverso do
quadrado da distância produz resultado que se assemelha ao resultado da
krigeagem proveniente de um modelo isotrópico. O ponto crítico, porém,
ocorre em regiões onde há superestimação e subestimação do teor de argila o
que difere do resultado oriundo da krigeagem, onde as mudanças de valores do
teor de argila são graduais.
· Os resultados produzidos pelos métodos de estimação média simples e vizinho
mais próximo (Figura 4.31) são menos expressivos com relação aos demais. O
método da média simples produz resultado que apresenta imbricação,
principalmente na região central da área de estudo (Figura 4.30). Já o método
de estimação relativo ao vizinho mais próximo (Figura 4.31), embora sendo o
que pior expressa a variabilidade espacial do fenômeno estudado, revela a
área de influência de cada ponto de observação. Tal informação é de grande
valia, como, por exemplo, numa análise preliminar para detecção de valores
amostrais suspeitos.
Uma forma alternativa de avaliar os resultados produzidos pelos respectivos
estimadores é através da operação de fatiamento aplicada às respectivas grades
geradas. Os resultados são apresentados a seguir.
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-35
De maneira análoga à análise anterior, constata-se que o resultado relativo ao fatiamento
da grade GKrV é o que melhor expressa a variabilidade do fenômeno, quando
comparado ao mapa de classes de solos concebido pelo CNPS-EMBRAPA (Figura
Geoestatística: Fundamentos e Aplicações
Geoprocessamento em Projetos Ambientais 5-36
4.34). Também verifica-se que o resultado referente ao fatiamento da grade GID (Figura
4.37) é próximo ao resultado do fatiamento da grade GKrI (Figura 4.36), mas apresenta
regiões que há superestimação e subestimação de valores. Os resultados referentes às
Figuras 4.38 e 4.39 são menos expressivos que os demais.
Este capítulo se encerra com a ilustração da Tabela 5.1, a qual apresenta, de
forma sintetizada, as principais razões pelas quais os resultados obtidos pela krigeagem
produziram resultados mais significativos em relação aos outros métodos de estimação.
TABELA 5.1 - CARACTERÍSTICAS DOS MÉTODOS DE ESTIMAÇÃO.
Métodos Geoestatísticos Métodos Convencionais
· os pesos são determinados a partir de
uma análise de correlação espacial
baseada no semivariograma.
li = f [g(h)]
· os pesos são determinados meramente
em função da distância.
li = f (di)
· Área de influência na interpolação é
indicada pelo alcance.
· raio de busca é arbitrário.
· Modela anisotropia, isto é, detecta as
direções de maior e menor
continuidade espacial do fenômeno.
· Anisotropia é ignorada.
· Trata redundância (“Clusters”), isto é,
atribui pesos adequados para
agrupamentos de amostras.
· Redundância é ignorada. Neste caso,
podem ocorrer superestimação ou
subestimação de valores.