SUMÁRIO
1 INTRODUÇÃO ______________________ __________________________________ _______________________ ______________________ ___________ 1 2 CONCEITOS BÁSICOS BÁSICOS _______________________ ___________________________________ _______________________ ______________ ___ 1 3 PREPARAÇÃO DE DADOS________________________ DADOS___________________________________ ______________________ ___________ 4 3.1 Composição de amostras de furos de sonda ______________________ _____________________________ _______ 4 3.1.1 Composição por bancadas _______________________ __________________________________ __________________ _______ 5
4 ANÁLISE ESTATÍSTICA ESTATÍSTICA _______________________ __________________________________ _______________________ ______________ 9 4.1 Conceitos de variáveis aleatórias e probabilidade______________________ probabilidade________________________ __ 10 4.2 Representações gráficas de variáveis aleatórias _______________________ _________________________ __ 12 4.3 Estatísticas descritivas de variáveis aleatórias ______________________ __________________________ ____ 15 4.3.1 Medidas de tendência central _________________ ____________________________ _____________________ __________ 15 4.3.2 Medidas de dispersão _______________________ __________________________________ _____________________ __________ 17 4.3.3 Medidas de forma ______________________ __________________________________ _______________________ _____________ __ 19 4.4 Modelos probabilísticos contínuos _______________________ __________________________________ _____________ __ 20 4.4.1 Distribuição normal ______________________________ _________________________________________ _______________ ____ 20 4.4.2 Distribuição lognormal ______________________ _________________________________ _____________________ __________ 22 4.5 Teorema do Limite Central ______________________ _________________________________ ___________________ ________ 23 4.6 Intervalo de confiança da média ________________________ ___________________________________ _____________ __ 25 4.7 Correlação e regressão ______________________ _________________________________ _______________________ ____________ 27
5 ANÁLISE GEOESTATÍSTICA GEOESTATÍSTICA _______________________ __________________________________ ___________________ ________ 29 5.1 Por quê variáveis regionalizadas ?___________________ ?______________________________ _________________ ______ 29 5.2 Variáveis regionalizadas _______________________ __________________________________ _____________________ __________ 31 5.3 O variograma ______________________ __________________________________ _______________________ ___________________ ________ 33 5.4 Relação entre semivariograma e a função covariância ___________________ ___________________ 35 5.5 Propriedades do variograma _______________________ __________________________________ _________________ ______ 37 5.6 Anisotropias _______________________ ___________________________________ _______________________ ___________________ ________ 38 5.7 Comportamento próximo à origem _______________________ ___________________________________ ____________ 39 5.8 Domínio do variograma _______________________ __________________________________ _____________________ __________ 40 5.9 Cálculo de variogramas experimentais ______________________ ________________________________ __________ 41 5.10 Modelos teóricos de variogramas___________ variogramas _______________________ _______________________ _____________ __ 42
6 ESTIMATIV ESTIMATIVAS AS POR KRIGAGE KRIGAGEM M ORDINÁRIA ORDINÁRIA ________ ____________ ________ ________ ________ _______ ___ 44 6.1 Definição da fronteira convexa__________________ convexa_____________________________ _____________________ __________ 44 6.2 Definição da vizinhança local _______________________ __________________________________ _________________ ______ 45
i
6.3 Definição da malha regular ______________________ _________________________________ ___________________ ________ 49 6.4 Krigagem ordinária___________ ordinária _______________________ ________________________ _______________________ _____________ __ 51 6.5 Validação cruzada ______________________ __________________________________ _______________________ _______________ ____ 64 6.6 Classificação de recursos/reservas minerais ______________________ ____________________________ ______ 67
7 ESTIMATIVAS POR COKRIGAGEM ORDINÁRIA____________ ORDINÁRIA _______________________ _____________ __ 71 7.1 Definições Básicas de Isotopia e Heterotopia ______________________ ____________________________ ______ 71 7.2 O variograma cruzado _______________________ __________________________________ _______________________ ____________ 71 7.3 O Modelo Linear de de Corregionalização ______________________ ________________________________ __________ 73 7.4 Cokrigagem ordinária _______________________ __________________________________ _______________________ ____________ 73
REFERÊNCIAS BIBLIOGRÁFICAS BIBLIOGRÁFICAS ______________________ _________________________________ _________________ ______ 76 ANEXO 1: DISTRIBUIÇÃO DISTRIBUIÇÃO NORMAL PARA X ENTRE 0 E 3,49 E AS INTEGRAIS Q(X) CORRESPONDENTES. CORRESPONDENTES. ______________________ _________________________________ _______________________ ____________81 81 ANEXO 2: VALORES CRÍTICOS CRÍTICOS DE T PARA ALGUNS ALGUNS NÍVEIS DE SIGNIFICÂNCIA ( IN KOCH & LINK, 1971 PÁG. 346)________________________ 346)________________________ 82
ii
1 INTRODUÇÃO
A geoestatística que foi definida inicialmente, por Matheron (1971, pág. 5), como a aplicação da Teoria das Variáveis Regionalizadas para a estimativa de depósitos minerais, tem hoje sua aplicação nas mais diversas áreas do conhecimento como: petróleo, hidrogeologia, meio ambiente, geotecnia, agronomia de precisão, oceanografia e reflorestamento. Como a geoestatística foi introduzida muito recentemente na grade curricular em cursos de graduação e de pós-graduação, há necessidade de promover cursos de extensão para disseminação da técnica, bem como para proporcionar uma reciclagem aos profissionais atuantes na área. Nesse sentido, surgiu a idéia de oferecer este curso, no qual introduziremos as técnicas e conceitos da geoestatística aplicada na análise e interpretação de dados, com o objetivo de fazer o melhor uso da informação disponível. Além disso, o planejamento deste curso, levou em consideração também à disponibilidade de um software totalmente nacional para que o aluno pudesse contar com uma licença acadêmica para que continuar seus estudos após o término do curso. Trata-se do sistema GeoVisual que foi desenvolvido para suportar o ensino de geoestatística em disciplinas de graduação e de pós-graduação ministradas regularmente no Instituto de Geociências – USP. 2 CONCEITOS BÁSICOS
A seguir vamos definir alguns conceitos básicos, cujo entendimento será de importância fundamental para interpretação dos resultados de uma análise geoestatística. Temos um problema a resolver, por exemplo, seja uma das seguintes questões: a) Qual é o teor médio de uma ocorrência mineral? b) Qual é o grau de contaminação por mercúrio no solo? c) Qual é a característica do solo para implantação de uma obra civil? A resposta para qualquer uma dessas questões está baseada na estimativa do atributo de interesse, através da amostragem. Am ostr os tr agem
Amostragem é o ato ou seleção de amostras como representativas do todo que se deseja estudar, tendo em vista sempre a limitação econômica do programa de amostragem. Normalmente, segundo Cochran (1963), as razões para a seleção de uma amostra são de ordem econômica - redução de custos - e apresentam como vantagens principais: maior rapidez, amplitude, flexibilidade e exatidão das 1
informações, em face da impossibilidade de se registrarem integralmente as especificações do mineral que se propõe conhecer. Base teórica
A teoria da amostragem é construída em torno do conceito que, se um número significativo de unidades representativas de uma população é selecionado sem enviesamento, o valor médio destas unidades irá aproximar-se da média da população (Barnes, 1980). Método Método s de amostragem
Existem basicamente três métodos de amostragem que poderiam ser considerados: a) aleatória simples; b) aleatória estratificada; c) sistemática. A amostragem sistemática, sempre que possível, é indicada para o cálculo de estimativas. Condição necessária
Estimativas de volumes, massas e teores devem ser baseadas em observações sistemáticas (amostragem sistemática) e interpretações da geologia (litologia e estrutura) e da mineralização (mineralogia, controles, distribuição e continuidades), segundo Vallée & Côte (1992). Fontes de erro erross
Os erros envolvidos na estimativa através da amostragem são devidos aos erros de amostragem e à variabilidade natural do fenômeno em estudo. Erros de amostr agem agem
“Qualquer amostragem - até mesmo a mais simples - comporta uma série de erros possíveis, alguns dos quais relacionados com a estrutura do minério, com sua distribuição e textura, outros decorrentes das técnicas usadas na amostragem, ou do modo como as técnicas são aplicadas, ou dos instrumentos de amostragem”, in Gy (1968). Este problema, infelizmente, não termina com a retirada da amostra, mas continua através da preparação, subdivisão e estágios de análise em laboratório, cada um dos quais é passível de erro, que podem 2
afetar a precisão ou influenciar sua exatidão. O erro de amostragem representa a diferença composicional entre a amostra de rocha e a parte do corpo rochoso com que se espera representá-la (Miesch, 1967). Segundo Koch & Link (1971), os erros de amostragem podem ser subdivididos, conforme a sua fonte de variação, em: a) erros de preparação: a preparação visa à redução do tamanho da amostra geológica, geralmente muito grande para fins de análise. Ela compreende uma série de operações não seletivas, tais como: redução da granulometria, mistura, homogeneização e subdivisão, cada uma das quais sujeita a erros; b) erros analíticos: são decorrentes da diferença entre o resultado da análise e a concentração na amostra original. Segundo Waeny (1979), os erros analíticos podem ser sistemáticos ou aleatórios. Os erros sistemáticos são aqueles que afetam as análises de maneira uniforme e decorrem da imperfeição dos instrumentos, da incorreção da técnica analítica, da impureza dos reagentes e de outros pequenos problemas, todos passíveis de controle ou atenuação (Waeny 1979). Os erros aleatórios, segundo o mesmo autor, são os de causa desconhecida, mas que podem ser localizados quando da análise periódica de um grupo de amostras. Independentemente do tipo de erro e do método analítico utilizado, existem limites de sensibilidade, além dos quais a determinação dos valores de concentração não é efetiva; c) erro total de amostragem: representa a soma dos erros decorrentes das etapas de amostragem e da preparação da amostra primária. Variabil idade natural
A variabilidade do fenômeno em estudo, que pode ser medida através do coeficiente de variação (razão da média pelo desvio padrão do conjunto de observações), pode representar a maior fonte de erro. O padrão de variação pode ser regular ou aleatório. Resultado
Os resultados da amostragem podem ser representados através de medidas diretas ou indiretas. As medidas diretas podem ser obtidas ‘in situ’ e/ou sobre a amostra e as medidas indiretas são obtidas através de sensores remotos (métodos geofísicos, sensoriamento remoto).
3
Natureza dos dados
Quanto à natureza os dados obtidos da amostragem podem ser classificados em qualitativos (cor, textura, etc.) e em quantitativos (composição química, teor, densidade, etc.). Inferência
A partir dos dados obtidos, os quais estão sujeitos a erros, devemos inferir as propriedades do todo amostrado ou do fenômeno em estudo. Se não houver variabilidade, poucas amostras serão suficientes para inferir o todo. Porém, na presença de variabilidade, a amostragem deve ser rigorosamente planejada para que, dentro da limitação econômica, seus resultados possam ser estendidos para o todo com um mínimo de erro. Aqui começa o problema para a geoestatística, pois a variabilidade entre as amostras deve ser reconhecida, mensurada e utilizada para posterior estimativa de porções não amostradas. Entretanto, antes de introduzir a geoestatística é necessário passar pela fase de preparação e analise estatística dos dados. 3 PREPARAÇÃO DE DADOS
As amostras podem ser coletadas de diversas formas (amostras de mão, fragmentos, canal, furos de sonda, etc.), nas quais deve-se garantir que apresentem o mesmo tamanho e massa. Por exemplo, amostras de mão de 2 kg, amostras de canal coletadas a cada 20 cm, testemunhos de sondagem a cada 2 m e, assim por diante. Porém, pelos motivos expostos a seguir, as amostras de furos de sonda necessitam de preparação para regularizar o intervalo de amostragem. 3.1 Composi ção de amostr as de furos de sonda
Geralmente o intervalo de amostragem nos furos de sonda não corresponde ao intervalo de trabalho na fase de avaliação de reservas, embora tenha sido necessário analisar as amostras segundo o intervalo de amostragem, sempre menor que o intervalo de trabalho. Justifica-se isto frente à necessidade de reconhecer e delimitar possíveis zonas ricas dentro da jazida. Além disso, as amostras individuais dos furos de sonda podem variar bastante em tamanho, comprimento e peso. Assim, a composição de amostras pelo agrupamento delas para o intervalo de trabalho, definido segundo a característica que se quer analisar, produzirá dados mais homogêneos e, portanto, com maior facilidade de interpretação. É importante especificar a unidade de amostragem utilizada na avaliação de reservas, pois, segundo Kim (1990), muitos problemas surgem pela não especificação da unidade de amostragem. Por exemplo, se uma jazida é avaliada com base na população de amostras de furos de sonda rotativa a 4
diamante, a produção da mina provavelmente não corresponderá às estimativas feitas, simplesmente porque a jazida não é lavrada com furos de sonda rotativa a diamante (Kim 1990). Segundo Barnes (1980), o objetivo de se fazer composições de amostras é obter amostras representativas de uma unidade mineralógica particular ou unidade de mineração. Segundo o mesmo autor, a unidade de amostragem é especificada no planejamento da amostragem e inclui o tamanho e modo de retirada física da amostra, assim como o intervalo de valores a ser analisado. O resultado da composição de amostras de furos de sonda é expresso como média ponderada do teor pelas espessuras selecionadas para o intervalo de trabalho, como mostra a equação a seguir: n
∑ t i ei
t c
= i =1n
(1)
∑ ei i =1
onde: n é o número de trechos para compor o intervalo de trabalho; t i é o teor do i-ésimo trecho; ei é a espessura do i-ésimo trecho. Existem muitos tipos de depósitos minerais, cada um dos quais irá requerer um tratamento especial dos dados amostrados para a obtenção dos melhores intervalos de composição para avaliação de depósito (Barnes, 1980). Basicamente, os tipos de composições possíveis em amostras de furos de sonda para o intervalo de trabalho são: -
bancadas; zona mineralizada.
Entretanto, vamos considerar apenas a composição por bancadas que é o caso mais comum de regularização de dados e não depende de interpretação prévia dos dados.
3.1.1 Composi ção por bancadas
O procedimento da composição por bancadas, é indicado para se fazer a avaliação de reservas em depósitos cuja lavra se dará a céu aberto. A composição por bancadas é feita aplicando-se a equação (1), onde as espessuras reais ou aparentes foram determinadas a partir de diferenças entre profundidades, dentro dos limites de cada bancada. No caso de furos inclinados, as espessuras são aparentes e o comprimento composto (CC) será maior que a altura da bancada, sendo tanto maior quanto menor a inclinação do furo, conforme ilustração na Figura 1. O comprimento composto (CC) pode ser calculado como: 5
CC =
altura da bancada
sen θ
Deve-se, em casos de furos inclinados, limitar a inclinação mínima que pode ser aceita para composições por bancadas, pois, por exemplo, para furos com inclinações de 30o, o comprimento composto será de 2 vezes a altura da bancada, como mostra a Tabela 1. Como recomendação, deve-se considerar esse ângulo mínimo igual a 20 o, pois este valor já dá um comprimento composto de aproximadamente três vezes a altura da bancada. A Figura 2 apresenta os fatores de multiplicação para inclinações variáveis entre 15 e 75 o.
Figura 1: Desenho esquemático mostrando o cálculo do comprimento composto em furos inclinados para cálculo de composições por bancada. Tabela 1: Fator de multiplicação da altura da bancada para cálculo do comprimento composto em furos inclinados. Inclinação do furo (o) 15 20 30 45 60 75
fator de multiplicação 3,864 2,924 2,000 1,414 1,547 1,035
Figura 2: Desenho ilustrando o problema da inclinação mínima de furos de sonda para a composição de amostras por bancadas. 6
Um outro problema relacionado à composição de furos inclinados está no cálculo das posições das novas amostras, ou seja, dos intervalos compostos. Deve-se estabelecer inicialmente se a posição é tomada no topo, meio ou pé da bancada. Recomenda-se utilizar sempre o meio da bancada como referência para localização das amostras compostas. Para o cálculo das posições das novas amostras deve-se determinar inicialmente o deslocamento horizontal (DH) das amostras (Figura 2-3), como a projeção do comprimento composto (CC), como segue: DH = CC cos(θ )
onde: θ é inclinação do furo de sonda em relação à horizontal. As posições das novas amostras nas bancadas podem ser calculadas recursivamente usando: X 1 Y 1
= X 0 + DH sen(φ )
= Y 0 + DH cos(φ )
onde: X1 e Y1 são as coordenadas da posição da nova amostra; X 0 e Y0 são as coordenadas da posição da amostra anterior e φ é o azimute do furo de sonda. A Figura 3 ilustra o procedimento do cálculo das posições das amostras compostas para a altura das bancadas.
Figura 3: Procedimento para cálculo das posições de amostras compostas para a altura das bancadas: representação do furo em seção (A) e projeção das coordenadas em planta (B). A Tabela 2 reproduz o log de um furo de sonda, com o qual pretende-se ilustrar o cálculo de composição por bancadas para o teor de Fe (%).
7
Tabela 2: Log do furo CP-23. Profundidade Descrição em metros De Até 0.00 8.35 Aterro da Estrada do Bota Fora 8.35 16.25 WH vermelho 16.25 28.00 WH amarelo com Goethita 28.00 32.43 HA pulverulenta 32.43 37.00 Rocha intrusiva 37.00 43.00 Filito amarelo 43.00 51.01 Filito cinza
Fe (%)
0.000 59.77 0 61.74 0 62.81 0 5.000 0.000 0.000
A Figura 4 apresenta esquematicamente o cálculo dos teores compostos de Fe para as bancadas 1410 e 1420 m.
Figura 4: Exemplo de composição de amostras por bancada para o furo CP-23 (vertical). Com o objetivo de exemplificar o cálculo de teores compostos por bancada para furos inclinados, considere-se o log de um furo de sonda, conforme os dados da Tabela 3. Tomando por base estes dados, o cálculo do teor composto de ferro para bancadas de 10 metros de altura será como exemplificado na Figura 5.
8
Tabela 3: Log do furo CP-62. Profundidade em metros de Até 0.00 4.20 4.20 16.40 16.40 26.40 26.40 28.70 28.70 29.72 29.72 33.04
Descrição
Fe (%)
Aterro 0.000 WH vermelho amarronzado 62.000 Hematita cinza escuro friável 60.000 Rocha intrusiva Hematita cinza escuro 60.000 Rocha intrusiva com Itabirito cinza 35.000 escuro 33.04 93.43 Itabirito cinza escuro, friável. 55.000
Figura 5: Exemplo de composição de amostras por bancada para o furo CP-62 (inclinado). 4 ANÁLISE ESTATÍSTICA
A análise estatística é uma etapa importante e deve preceder a análise geoestatística e a estimativa por krigagem ordinária. Esta análise permite sumariar os dados obtidos, conferir a base de dados e, inclusive, reconhecer valores anômalos. Trata-se em conhecer melhor os dados que serão utilizados para estimativas e inferências e, portanto, eles devem estar isentos de erros de digitação. Antes de introduzir os conceitos estatísticos, seria interessante rever brevemente alguns conceitos de variáveis aleatórias e probabilidade. 9
4.1 Conceitos de variáveis aleatórias e prob abilid ade
Vamos introduzir e rever alguns conceitos importantes sobre variáveis aleatórias e probabilidade, antes de passar à análise estatística propriamente dita. Variáveis aleatórias
Uma variável cujo valor é determinado pela realização de um experimento é denominada variável aleatória. As variáveis aleatórias podem ser subdivididas em duas classes: discretas e contínuas. Variáveis aleatórias di scretas
Uma variável aleatória discreta é aquela que tem um número contável de realizações possíveis. Exemplo: cor do solo, segundo uma escala padrão de cores. Variáveis aleatórias con tínuas
Uma variável aleatória contínua é aquela que pode assumir qualquer valor num segmento contínuo da linha dos números reais. Exemplo: teores, massas, volumes são geralmente medidos em uma escala contínua. Notação
Letras maiúsculas serão usadas para referir as variáveis aleatórias. Exemplo: X, representando o teor de sílica em minério de ferro. Letras minúsculas referem-se a um valor específico da variável aleatória. Exemplo: x=3.23 %, é o teor de sílica em uma amostra específica. Probabili dade de realizações igualmente possíveis
Se A e B são eventos do espaço amostral S, então a probabilidade do evento A, P(A) é igual à freqüência do evento A (realizações simples igualmente possíveis) sobre n o tamanho da amostra: P( A) =
f A n
10
e, similarmente: P( B ) =
f B n
Então: 1) 0 ≤ P( A) ≤ 1 e 0 ≤ P(B ) ≤ 1 2) P(S ) = 1 3) Se A e B são eventos mutuamente exclusivos, ou seja, não podem ocorrer simultaneamente (A∩B=0), então P( A ∪ B ) = P( A) + P( B ) . Distribuição de probabilidade
As probabilidades p(x1), p(x2), ... ,p(xn) associadas aos valores possíveis x1, x 2, ... ,xn de uma variável aleatória X constituem a distribuição de probabilidade de X.
Função de probabilidade
Ao conjunto de pares ordenados (x i, p(xi)) denomina-se função de probabilidade de X. Função densidade de probabilidade
Determina as probabilidades teóricas associadas às variáveis aleatórias contínuas. A função densidade de probabilidade permite calcular a probabilidade da variável aleatória estar no intervalo finito [a,b] como segue: b
P(a ≤ x ≤ b ) = ∫ f ( x')dx' a
Função de di stribui ção acumulada
Dá a probabilidade para a variável aleatória x’ ser menor ou igual a x: P( x' ≤ x ) = F ( x ) =
x
∫ f ( x')dx'
−∞
11
4.2 Representações gráficas de variáveis aleatórias
As observações amostradas de uma variável aleatória podem ser representadas graficamente com o objetivo de estudar a sua distribuição dentro do intervalo amostrado. As representações gráficas mais utilizadas são o histograma e a curva acumulativa. Histograma
A análise estatística começa pelo estudo da distribuição de freqüências, a qual descreve como as unidades de uma população estão distribuídas sobre o intervalo amostrado. A distribuição de freqüências pode ser do tipo simples ou acumulada. A distribuição de freqüências do tipo simples é construída tabulando-se os dados de alguma característica medida do depósito (teor, espessura, etc.) em intervalos constantes; os dados assim agrupados podem ser representados graficamente na forma de histograma, lançando-se os intervalos de medida em abscissa e as freqüências em ordenada. A Figura 6 apresenta um histograma para uma variável aleatória tipicamente encontrada na análise de dados geológicos. O histograma proporciona uma representação gráfica que permite visualizar a distribuição dos dados. 25
%
20 15 10 5 0 0.5
1.0
1.5
2.0
Valores dos dados
Figura 6: Histograma para uma variável aleatória típica de dados geológicos.
Curva acumulativa
O procedimento para obtenção de freqüências acumuladas é o mesmo que o do tipo anterior, porém as freqüências dos dados agrupados nos intervalos são agora acumuladas. A curva acumulativa é a representação gráfica obtida do lançamento das freqüências acumuladas em ordenada e os intervalos de medida 12
em abscissa. A simples união dos pontos sobre a curva acumulativa, com segmentos de reta, gera o polígono de freqüências acumuladas. Entretanto, este tipo de representação não tem sido utilizado, pois não permite determinar com precisão o valor da variável de interesse para um determinado percentil, ou viceversa. Assim, mais recentemente as curvas acumulativas têm sido construídas a partir das freqüências acumuladas de todos os dados do conjunto amostrado. Nesse caso, as freqüências simples são calculadas para todos os valores observados individualmente, atribuindo-se uma freqüência igual a 1/n. As freqüências simples assim obtidas podem então ser acumuladas gerando as freqüências acumuladas. As freqüências acumuladas são geralmente lançadas em escala de probabilidade aritmética, que tem a propriedade de identificar graficamente se a distribuição em estudo segue uma distribuição normal (os pontos deverão estar alinhados sobre uma reta). A escala de probabilidade aritmética é utilizada freqüentemente para a representação gráfica de distribuições de freqüências acumuladas, pois permite determinar rapidamente se a distribuição em estudo é normal ou lognormal. Esta escala é obtida da integração da função densidade de probabilidade da distribuição normal de - ∞ a X: P( X ) =
X
∫ f ( x)dx
−∞
P(X) é a probabilidade acumulada de - ∞ a X, ou a uma porcentagem em relação à área total da curva, e corresponde à porcentagem acumulada na escala de probabilidade aritmética, como ilustrado na Figura 7. Se a distribuição de freqüências acumuladas for do tipo normal ou aproximar-se dele, os pontos deverão alinhar-se numa reta, quando lançados em gráfico de probabilidade aritmética (abscissa em escala aritmética). Contudo, se os pontos desenharem um “S”, significa que a distribuição não é normal e deve ser testada a hipótese de lognormalidade, lançando-se os pontos em gráfico de logprobabilidade aritmética (abscissa em escala logarítmica). Se os pontos alinharem-se em torno de uma reta significa que a distribuição é lognormal. Assim, é possível verificar graficamente se a distribuição é normal ou lognormal, representando a distribuição de freqüências acumuladas em gráficos de probabilidade ou logprobabilidade aritmética, respectivamente.
13
Figura 7: A escala de probabilidade aritmética resultando da integração da função densidade de probabilidade da distribuição normal. A Figura 8 apresenta uma curva acumulativa para a mesma variável aleatória representada no histograma da Figura 6. 99.99
A D A 99.95 L 99.90 U M U99.50 C99.00 A %
95.00 90.00 80.00 70.00 60.00 50.00 40.00 30.00 20.00 10.00 5.00 1.00 0.50 0.10 0.05 0.01 0.5
1.0
1.5
2.0
Valores dos dados
Figura 8: Curva acumulativa para os valores da variável aleatória típica representada no histograma da Figura 6.
14
4.3 Estatísticas descrit ivas de variáveis aleatórias
As distribuições de freqüência apresentam características intrínsecas ao fenômeno em estudo, por exemplo, com menor ou maior variabilidade. As estatísticas descritivas são utilizadas para caracterizar numericamente as distribuições de freqüência. Estas estatísticas podem ser obtidas através de: medidas de tendência central, medidas de dispersão e medidas de forma da curva.
4.3.1 Medidas de tendência central
Ao se estudar uma distribuição de freqüências, o objetivo é determinar o valor mais provável dessa distribuição, que pode ser encontrado a partir de medidas de tendência central: média, mediana e moda. Média
A média ou esperança matemática de uma variável aleatória X é definida como: n
E [ X ] = ∑ xi p( xi ) i =1
onde: p( xi ) é a probabilidade associada à ocorrência de xi . Se os valores x1 , x2 ,L, xn representam os valores possíveis e estes são igualmente possíveis ( p( x1 ) = p( x2 ) = L = p( xn )) , então a média torna-se: n
1
i =1
n
E [ X ] = X = ∑ xi
=
1
n
∑ xi
n i =1
Propri edades da média
Tem-se a seguir algumas propriedades associadas à média (Fonseca & Martins, 1982, págs. 40-41): a) a média de uma constante é a própria constante; 15
E [K ]
=
K
b) a média de uma variável aleatória X, multiplicada por uma constante, é igual a constante multiplicada pela média de X; E [KX ]
=
K . E [ X ]
c) a média da soma ou diferença de duas variáveis aleatórias é a soma ou diferença das médias; E [ X ± Y ]
= E [ X ] ± E [Y ]
d) a média de uma variável aleatória somada ou subtraída de uma constante é igual à média dessa variável somada ou subtraída da mesma constante; E [ X ± K ]
= E [ X ] ± K
Observação: esta propriedade é particularmente importante para transformação de variáveis visando ajustar a sua média sem, contudo, alterar a variância. e) a média de uma variável aleatória subtraída de sua própria média é zero; E X − X
=0
Mediana
A mediana corresponde ao valor da variável aleatória a 50% da distribuição acumulada de freqüências. No exemplo da Figura 7, o valor da variável aleatória correspondente a 50% é igual a 1,204. Moda
A moda corresponde à classe mais freqüente verificada no histograma. Para o exemplo da Figura 6, a moda corresponde à classe 1,2 – 1,3. Quando duas modas são verificadas, com as mesmas freqüências, então a distribuição é bimodal.
16
4.3.2 Medidas de dispersão
Da mesma forma que existem várias maneiras para medir a tendência central dos dados, há também várias maneiras para medir a dispersão em torno da média: variância e desvio padrão, coeficiente de variação e Teorema de Chebyshev.
Variância e desvio padrão
A dispersão dos valores em torno da média é medida pela variância, que é determinada como: Var [ X ]
n
= ∑ ( xi − X ) 2 p( xi ) i =1
Novamente, assumindo as probabilidades de ocorrência dos n valores possíveis iguais entre si, ou seja, iguais a 1/n, tem-se: Var [ X ]
=
2
S
n
1
i =1
n
= ∑ ( xi − X ) 2
=
1
n
∑ ( xi − X ) n
2
i =1
que é a equação usual da variância, ou desenvolvendo-a, tem-se: 2
S
= E [ X 2 ] − E [ X ]2
O desvio padrão é simplesmente a raiz quadrada da variância e é expresso na mesma unidade dos valores originais. Propriedades da variância
As propriedades associadas à variância, segundo Fonseca & Martins (1982) são: a) a variância de uma constante é zero; Var [K ]
=
0
17
b) a variância de uma variável aleatória multiplicada por uma constante é igual a variância da variável aleatória multiplicada pelo quadrado dessa constante; Var [KX ]
=
K Var [ X ] 2
c) a variância de uma variável aleatória somada ou subtraída de uma constante é igual à variância da variável aleatória; Var [ X ± K ] = Var [ X ]
d) a variância da soma ou diferença entre duas variáveis aleatórias independentes é a soma das respectivas variâncias. Var [ X ± Y ]
= Var [ X ] + Var [Y ]
Coeficiente de vari ação
O coeficiente de variação, que é uma outra medida de dispersão, é obtido pela divisão do desvio padrão pela média: CV =
S X
Como o coeficiente de variação é adimensional, ele é freqüentemente utilizado para comparar a dispersão relativa de valores em torno da média entre diferentes distribuições, como, por exemplo, para comparação e classificação de depósitos minerais segundo a variabilidade natural, medida por meio desta estatística. Teorema de Chebyshev
De acordo com o Teorema de Chebyshev, para qualquer função densidade de probabilidade, a proporção da variável aleatória dentro de ± K desvios padrão em torno da média é sempre no mínimo 1 −
1 K 2
, onde K é qualquer número
positivo maior que 1. A Tabela 4 mostra a proporção de Chebyshev para alguns valores de K.
18
Tabela 4: Proporção de Chebyshev da variável aleatória estar dentro de ± K desvios padrão em torno da média. K
2 3 4
≥ 1−
1 K 2
0.75 0.89 0.94
4.3.3 Medidas de forma
As distribuições de freqüências podem ser caracterizadas também quanto à forma, através das medidas de assimetria e curtose. As simetria
Assimetria é a medida do grau de simetria de uma distribuição de freqüências em torno da média, a qual pode apresentar uma assimetria positiva se a cauda da distribuição estiver à direita da média e negativa se estiver à esquerda. A assimetria positiva é observada na maioria das distribuições de freqüências de variáveis de depósitos minerais com alta variabilidade natural (metais raros, ouro, urânio, etc.). O coeficiente de assimetria pode ser calculado a partir do terceiro momento centrado na média: n
CA=∑ ( xi
− X ) 3 / S 3
i =1
Curtose
A curtose é a medida do grau de achatamento de uma distribuição em relação à distribuição normal (Spiegel, 1967), que reflete a dispersão dos valores em torno da média. O coeficiente de curtose é calculado a partir do quarto momento em torno da média: n
CC =∑ ( xi
− X ) 4 / S 4
i =1
19
4.4 Modelos probabilísticos contínuos
Serão apresentados neste item os principais modelos probabilísticos utilizados para descrever o comportamento de variáveis aleatórias contínuas encontradas na análise de dados geológicos. Tais modelos são representados pela distribuição normal e lognormal.
4.4.1 Distribuição normal
A distribuição normal ou gaussiana é a mais comumente utilizada em estatística, pois sob esta forma de distribuição de freqüências encontra-se um grande número de variáveis aleatórias em muitos campos de aplicação. A função densidade de probabilidade, que descreve matematicamente esta distribuição, é dada por:
f ( x)=
1 σ 2π
1 / 2[( x − µ ) / σ ] 2
e−
(2)
onde: f ( x ) é a função densidade de probabilidade; x é uma observação; µ e σ são respectivamente a média e o desvio padrão que definem a forma da curva. A Figura 9 apresenta os gráficos da função densidade de probabilidade, nos quais encontram-se delimitadas as áreas correspondentes a 68, 95 e 99,7% da distribuição, ou seja, equivalentes aos intervalos µ±σ, µ±2σ e µ±3σ, respectivamente.
Figura 9: Gráficos da função densidade de probabilidade da distribuição normal para áreas correspondentes a 68% (A); 95% (B) e 99,7%(C) de distribuição. A Tabela 5 compara as proporções de Chebyshev para qualquer variável aleatória com a variável aleatória normal.
20
Tabela 5: Valores de probabilidade que x está no intervalo ± K desvios padrão em torno da média (Exell, 1998). K
1 2 3 4
Qualquer variável aleatória n.d. ≥ 0.75 ≥ 0.8889 ≥ 0.9375
Variável aleatória normal 0.68 0.95 0.997 0.99994
Desta tabela pode-se concluir que de todas as variáveis aleatórias possíveis com a mesma variância, a variável aleatória contínua é a que está mais concentrada em torno da média (Exell, 1998), daí a sua grande utilidade na prática. As áreas sob a distribuição normal podem ser facilmente calculadas integrando-se a função densidade de probabilidade [equação (2)], por exemplo, da posição X até +∞, resultando na área Q(X), como ilustrado na Figura 10. Assim, pode-se calcular as áreas sob a distribuição normal entre 0 e 3,49 resultando numa forma de tabela da distribuição normal (Anexo 1). A distribuição normal é, sem dúvida, a distribuição teórica mais utilizada na prática, pois é matematicamente conveniente de se trabalhar com ela, uma vez que suas propriedades são bastante conhecidas. Em geral, a grande maioria das variáveis aleatórias segue uma distribuição normal ou, no mínimo aproximadamente normal. Mesmo para observações que não apresentam uma distribuição normal, esta pode ser aplicada se o problema puder ser resolvido considerando o comportamento de uma variável formada pelo cálculo da estatística de um conjunto de observações, ao invés das observações originais (vide Teorema do Limite Central).
Figura 10: Gráfico da distribuição normal mostrando a área correspondente a integral da função densidade de probabilidade de X à +∞.
21
4.4.2 Distribuição lognormal
A distribuição lognormal é um tipo encontrado em muitos problemas de avaliação de reservas (principalmente em casos de metais raros), caracterizandose por uma distribuição com assimetria positiva, onde ocorre uma grande quantidade de valores baixos e uns poucos valores altos. Formalmente, a distribuição lognormal é definida como uma distribuição contínua caracterizada pela propriedade dos logaritmos das observações seguirem uma distribuição normal (Koch & Link, 1971). A função densidade de probabilidade da distribuição lognormal é dada por: 1 1 / 2[(log x −α ) / β ] 2 f ( x)= e− x β 2π
(3)
onde: α é a média dos logaritmos de x; ß é o desvio padrão dos logaritmos de x em relação a α. Os dois parâmetros α e ß2 definem a forma da curva de distribuição de probabilidades. A distribuição lognormal é sempre assimétrica para a direita (assimetria positiva), sendo que o grau de assimetria depende somente do valor de ß2, que corresponde à variância dos logaritmos das observações (Koch & Link, 1971). O objetivo básico da transformação não linear (logarítmica) observada na equação (3) é, segundo Koch & Link (1971), a mudança da forma da distribuição para uma distribuição normal ou aproximadamente normal. Entretanto, algumas distribuições de freqüências apresentam-se após a transformação logarítmica com certa assimetria negativa (assimétrica para a esquerda), que pode ser corrigida pela adição de uma constante - o terceiro parâmetro - às observações originais, antes da transformação logarítmica (Krige, 1978). A função densidade de probabilidade da distribuição lognormal a três parâmetros é descrita por: f ( x)=
1 ( x + C ) β 2π
1 / 2[(log( x +C ) −α ) / β ]2
e−
A constante C - terceiro parâmetro da distribuição lognormal - pode ser estimada, segundo Landim (1985), como: C =
− p1 p 2 p1+ p 2−2 M M
2
22
onde: M é a mediana, ou seja, o valor de teor correspondente a 50% da distribuição; p1 é o valor de teor correspondente a um percentil entre 5 e 20%; p2 é o valor de teor correspondente a um percentil entre 100-p1. Na Figura 11, tem-se as curvas de distribuição lognormal com α e C iguais a zero e para três valores de ß 2 (2, 0,5 e 0,1).
Figura 11: Curvas de distribuição lognormal com α e C iguais a zero e três valores de β2, segundo Aitchison & Brown (1957, apud Koch & Link, 1971).
4.5 Teorema do Limite Central
Segundo Barnes (1980), o Teorema do Limite Central é um dos mais importantes teoremas da estatística matemática relacionada a distribuições de freqüências de amostragem e pode ser enunciado como: "Se amostras aleatórias de tamanho fixo são retiradas de uma população cuja distribuição teórica é de forma arbitrária, mas com média e variância finitas, a distribuição das amostras tende mais e mais a uma distribuição normal com média µ e variância σ2/n tanto quanto o tamanho das amostras aumenta”. Se X1, X2, ..., X n são valores de uma variável aleatória, a média X é: X =
( X 1 + X 2 + ... + X n ) n
O valor esperado de X , usando as propriedades (b) e (c) da média, é:
[ ] = 1n ( E [ X ] + E [ X ]+ .... + E [ X ])
E X
1
2
n
ou 23
[ ]=
E X
1 n
( n.µ )
=
µ
onde µ é a média populacional A variância de X é calculada como: Var X
= Var [( X 1 + X 2 + .... + X n ) / n]
aplicando-se as propriedades (b) e (d) da variância tem-se: 2
⎛ 1 ⎞ Var [ X ] = ⎜ ⎟ (Var [ X 1 ] + Var [ X 2 ] + ....Var [ X n ]) ⎝ n ⎠ ou
[ ] = ⎛ ⎜ 1 ⎞⎟nσ
Var X
2
⎝ n ⎠
2
= σ 2 / n
A distribuição de X tem média µ e variância σ2/n que se aproxima da distribuição normal tanto quanto aumenta o tamanho da amostra. Na maioria dos casos a aproximação é boa a partir de 40 amostras. A Figura 12 ilustra muito bem o que enuncia o Teorema do Limite Central. Amostras aleatórias de tamanho fixo são retiradas de distribuições arbitrárias (Figura 12A), a partir das quais tem-se: as distribuições das médias para 2 amostras (Figura 12B), para 4 amostras (Figura 12C) e para 25 amostras (Figura 12D). Observe-se que a média das amostras tende a µ e variância σ2/n, tanto quanto aumenta o tamanho das amostras. As médias praticamente permanecem, enquanto as variâncias diminuem na proporção da raiz quadrada do número de amostras. As n amostras da variável aleatória X são, na verdade, em problemas de avaliação de reservas, os teores de n blocos de cubagem que compõem o depósito mineral em avaliação. O teor médio no bloco de cubagem é determinado como a média ponderada dos teores de amostras de furos vizinhos. O ponderador variará de acordo com o método escolhido na avaliação de reservas, seja ele convencional ou computacional. Portanto, o teor médio do depósito será igual à média dos teores calculados nos blocos de cubagem, como está assegurado pelo Teorema do Limite Central. Por isso, o teor médio do depósito, desde que as informações coletadas estejam bem representadas no mesmo, pode ser determinado com precisão, pois mesmo com a alteração da dimensão do bloco, o teor médio deverá se manter na mesma faixa de valores, enquanto a variância diminuirá com o tamanho do bloco de cubagem.
24
4.6 Intervalo de conf iança da média
Calculado o valor médio da variável de interesse, pode-se determinar o intervalo de confiança associado ao mesmo, a um determinado nível de confiança. O intervalo de confiança pode ser calculado a partir da estatística t: t =
X − µ S n
(A)
(B)
(C)
(D)
Figura 12: Populações arbitrárias (A), das quais são selecionadas aleatoriamente amostras de tamanho fixo: n=2 (B); n=4 (C) e n=25 (D), segundo Lapin (1982, apud Davis, 1986). 25
que tem uma nova distribuição de amostragem. A expressão para o intervalo de confiança da média populacional ao nível de confiança de 90% é: X − t 5%
s n
< µ < X + t 5%
s n
A distribuição t é simétrica e depende somente do número de graus de liberdade, que no caso, da estimativa da média, é igual a n-1 (n=número de amostras). O Anexo 2 apresenta os valores críticos de t para alguns níveis de significância. Quando o número de graus de liberdade tende ao infinito, a distribuição t tende à distribuição normal, como está ilustrado na Figura 13. Observe-se que o valor crítico de t com graus de liberdade tendendo ao infinito, a um nível de significância de 10%, corresponde à variável aleatória padronizada da distribuição normal para uma área equivalente a 10% (Anexo 1).
Figura 13: Distribuição t de Student para vários graus de liberdade. 26
4.7 Correlação e regressão
Muitas vezes é necessário estudar a relação mútua entre duas variáveis aleatórias (X e Y) com o objetivo de verificar se elas encontram-se correlacionadas ou não, por exemplo, para o cálculo de co-estimativas. Para isso necessita-se das estatísticas que correlacionem duas variáveis aleatórias. Covariância
A variância mede a dispersão de uma variável aleatória X em torno da sua média X , enquanto a covariância mede a dispersão ou como encontram-se correlacionadas duas variáveis aleatórias simultaneamente: Cov( X , Y ) = E [( X − X )(Y − Y )] ,
que pode ser desenvolvida como: Cov( X , Y ) = E [ XY ] − E [ X ] E [Y ]
Observe-se que ao contrário da variância que é sempre positiva, a covariância pode ser positiva ou negativa, dependendo da relação existente entre as variáveis aleatórias. Coeficiente de corr elação
O resultado da covariância nem sempre é de fácil interpretação, pois depende dos valores associados às variáveis aleatórias X e Y. Assim, comumente utiliza-se de outra medida derivada da covariância denominado coeficiente de correlação, conforme a seguinte expressão: Corr ( X , Y ) =
Cov( X , Y ) Var ( X )Var (Y )
,
que tem a vantagem de estar normalizado no intervalo –1 a +1. Um valor próximo de zero indica a falta de correlação entre as duas variáveis.
27
Reta de regressão
Dado um conjunto de pares ordenados {x i, y i, i=1,n}, pode-se determinar a relação funcional y = f ( x ) , através do método dos mínimos quadrados. Vamos exemplificar neste item a obtenção da reta dos mínimos quadrados: y ∗ = a + bx . Segundo este método, a diferença elevada ao quadrado entre os valores observados (y) e calculados (y*) deve ser a mínima possível. Assim, devemos minimizar: n
(
− yi∗ ) , 2
S = ∑ yi i =1
onde yi∗ = a + bxi e, portanto: n
S = ∑ ( yi
− a − bxi )2
i =1
Para encontrarmos o mínimo de S com relação aos coeficientes (a, b), calculamos as derivadas parciais e igualamos a zero: dS da dS db
n
= ∑ 2( yi − a − bxi )(− 1) = 0 i =1
n
= ∑ 2( yi − a − bxi )(− xi ) = 0 i =1
desenvolvendo, tem-se: n
n
i =1
i =1
n
a ∑ 1 + b∑ xi
n
a ∑ xi i =1
= ∑ yi i =1
n
n
i =1
i =1
+ b∑ xi2 = ∑ xi yi
Por fim, os coeficientes procurados são: n
∑ yi
a
= i =1
b=
n
− b∑ xi i =1
n
Cov( X , Y ) Var [ X ]
28
5 ANÁLISE GEOESTATÍSTICA
Na estatística trabalhamos com realizações de variáveis aleatórias; na geoestatística trabalhamos com as funções aleatórias onde as amostras são vistas como realizações de uma variável aleatória que, por sua vez, é função das coordenadas espaciais. A geoestatística envolve a análise e predição de fenômenos espaciais ou temporais, tais como: teores de minério, porosidades, concentração de poluentes, preço do petróleo no tempo, etc. À etapa de estudo e modelagem da correlação espacial denomina-se análise geoestatística. É desta análise que se obtém a ferramenta básica da estimativa por meio da krigagem ordinária, que é o variograma. Após a análise geoestatística, pode-se fazer predições ou simulações estocásticas em pontos não amostrados para melhor compreensão do fenômeno espacial em estudo. Com o modelo de variograma reconhecem-se anisotropias (feição particular dos métodos geoestatísticos), bem como uma idéia da variabilidade a pequenas distâncias dada pelo comportamento próximo à origem. Cabe salientar que a krigagem, como método de estimativa da variável de interesse, só deve ser utilizada quando o variograma experimental for estruturado, ou seja, se a variabilidade não for totalmente aleatória (efeito pepita puro).
5.1 Por quê variáveis regionalizadas?
As variáveis regionalizadas, que representam os valores de variáveis referenciadas geograficamente, foram introduzidas para descrever quantitativamente variações espaciais em corpos de minério. Este item, baseado no trabalho de Royle (1979), justifica porque as variáveis regionalizadas são dependentes de suas posições espaciais relativas e também mostra como se pode medir as variações espaciais. Para o desenvolvimento deste item, serão consideradas as seguintes séries de números: Série A:
1
7
3
6
2
9
4
8
5
Série B:
1
3
5
7
9
8
6
4
2
As características estatísticas dessas duas séries de números, medidas através da média e variância (Tabela 6), são idênticas, pois apresentam os mesmos valores. Entretanto, essas duas séries são bem diferentes, pois resultam de dois tipos distintos de mineralização. Tabela 6: Estatísticas medidas para as duas séries de números. Série A B
Média 5 5
Variância 6,67 6,67 29
Assim, as estatísticas, obtidas por métodos clássicos, não conseguem reconhecer a diferença existente entre as duas séries em estudo, pois consideram as amostras independentes entre si. Por outro lado, se fosse considerada a posição espacial relativa de cada amostra, poder-se-ia distinguir as duas séries de números. Uma possibilidade seria medir a diferença entre os valores de amostras separadas por uma determinada distância. Como a simples soma das diferenças tenderia a anular-se, optou-se pela soma do quadrado das diferenças, que dividido pelo número de pares dá sentido a uma medida de variância, com significado espacial, pois é dependente da distância utilizada. A variância espacial pode ser calculada para várias distâncias, ou para vários intervalos de amostragem, como segue: variância espacial para um intervalo de amostragem: 2
+ (7 − 3)2 + (3 − 6)2 + (6 − 2)2 + (2 − 9)2 + (9 − 4)2 + (4 − 8)2 + (8 − 5)2
/ 8 = 22
2
+ (3 − 5)2 + (5 − 7)2 + (7 − 9)2 + (9 − 8)2 + (8 − 6)2 + (6 − 4)2 + (4 − 2)2
/ 8 = 3,63
A : (1 − 7 ) B : (1 − 3)
variância espacial para dois intervalos de amostragem: A : (1 − 3)
2
+ (3 − 2)2 + (2 − 4)2 + (4 − 5)2 + (7 − 6)2 + (6 − 9)2 + (9 − 8)2
/7 = 3
2
+ (5 − 9)2 + (9 − 6)2 + (6 − 2)2 + (3 − 7 )2 + (7 − 8)2 + (8 − 4)2
/ 7 = 12,86
B : (1 − 5)
e, assim sucessivamente. Calculando-se a variância espacial até quatro intervalos de amostragem tem-se os resultados mostrados na Tabela 7. Tabela 7: Variâncias espaciais para as séries A e B, determinadas até quatro intervalos de amostragem. intervalo de amostragem 1 2 3 4
var. esp. A
var. esp. B
22,00 3,00 23,67 3,80
3,63 12,86 23,83 29,60
Os dados da Tabela 7, podem ser representados sob forma gráfica, lançando-se as variâncias espaciais em função dos intervalos de amostragem, como está mostrado na Figura 14.
30
Figura 14: Variância espacial em função dos intervalos de amostragem para as séries A e B. Observando-se o gráfico da Figura 14, verifica-se que a série A é muito errática, enquanto a série B é mais uniforme. Na série B, as variâncias espaciais aumentam conforme o intervalo de amostragem, pois a correlação entre os valores diminui com a distância. Esse comportamento seria desejado em todos os corpos de minério, dentro de alguma escala de amostragem. Uma maneira prática para verificar se há correlação espacial nos dados, segundo Bon (1979), é calcular a variância espacial para um intervalo de amostragem e comparar com a variância amostral; se a variância espacial for menor, então há correlação, caso contrário não há. Veja, por exemplo, que na série B, a variância espacial para um intervalo de amostragem é igual a 3,63 para uma variância amostral de 6,67, portanto, com boa correlação espacial, enquanto para a série A, a mesma variância é igual a 22,00, ou seja, não apresenta correlação espacial.
5.2 Variáveis regionalizadas
Uma variável regionalizada é qualquer função numérica com uma distribuição espacial, que varia de um lugar a outro com continuidade aparente, mas cujas variações não podem ser representadas por uma função determinística (Blais & Carlier, 1968 apud Olea, 1975). O termo variável regionalizada foi escolhido por Matheron (1965, apud Huijbregts, 1975) para enfatizar as feições particulares dessas variáveis. Em geologia, todas as observações quantitativas feitas em duas ou três dimensões (área ou volume, respectivamente), sejam elas geoquímicas, geofísicas, sedimentológicas, etc., podem ser consideradas como exemplos de variáveis regionalizadas. A definição de uma variável regionalizada como uma variável distribuída no espaço é puramente descritiva e não envolve qualquer interpretação probabilística. Uma variável aleatória é aquela que recebe um certo número de valores, de acordo com uma certa distribuição de probabilidades (Journel & Huijbregts, 1978). O teor de um elemento num ponto x 1 do depósito pode ser considerado como uma realização particular de uma variável aleatória Z(x 1) 31
definida no ponto x1. Segundo Journel & Huijbregts (1978), denomina-se função aleatória Z(x) o conjunto de teores Z(x) para todos os pontos x dentro do depósito [i.e. variável regionalizada Z(x)]. A interpretação probabilística de uma variável regionalizada, como uma realização particular de uma certa função aleatória Z(x), tem um significado operacional quando for possível inferir toda ou parte da lei de probabilidades que define essa função aleatória na sua totalidade (Journel & Huijbregts, 1978). A maioria das variáveis regionalizadas apresenta um aspecto aleatório, consistindo de variações altamente irregulares e imprevisíveis, e um aspecto estruturado, refletindo as características estruturais do fenômeno regionalizado (Kim, 1990). Uma formulação apropriada para solução de problemas de estimativa deve levar em consideração essas duas características aparentemente contraditórias, por meio de uma representação simples da variabilidade espacial (Journel & Huijbregts 1978). A Teoria das Variáveis Regionalizadas tem por objetivos o estudo e representação das propriedades estruturais das variáveis regionalizadas para resolução de problemas de estimativa. Hipótese intrínseca
"Um conceito básico na Teoria das Variáveis Regionalizadas é a chamada hipótese intrínseca, a qual implica que uma função (a função intrínseca) descreve o comportamento espacial da variável regionalizada dentro do espaço e que essa função é uma característica intrínseca da regionalização. A função intrínseca é na verdade o chamado semivariograma. Em outras palavras, a geoestatística assume que a distribuição das diferenças entre dois pontos amostrais (estatística de dois pontos) é a mesma para todo o depósito e que ela depende apenas da distância e orientação entre os pontos. Essa é a conceituação geoestatística da hipótese intrínseca, algumas vezes referenciada como hipótese de quaseestacionaridade. A variação espacial é estacionária se ela puder ser reconhecida em todas as partes do espaço, ou seja, o variograma é o mesmo onde quer que se amostre. A estacionaridade usada na Teoria das Variáveis Regionalizadas é a estacionaridade de segunda ordem das diferenças entre a variável Z(x) e a variável Z(x+h) nos pontos (x) e (x+h), onde (Z(x); x ∈ D) é um processo estocástico a valores reais, definido sobre um domínio D em R, R2 ou R3", in IPT (1989). Características qualitativas das variáveis regionalizadas
Segundo Bubenicek & Haas (1969), as características qualitativas de variáveis regionalizadas, que os métodos estatísticos convencionais não conseguem reconhecer, são:
32
Localização
Os valores de uma variável regionalizada são dependentes de suas funções espaciais relativas dentro do campo geométrico (depósito). Além disso, estes valores são dependentes do tamanho da amostra, forma e orientação (suporte amostral); Suporte
Por vezes a variável regionalizada Z(x) não está definida num ponto, mas sobre uma área ou volume centrado em x. A unidade amostral básica sobre a qual a variável é medida chama-se suporte (IPT, 1989); Continuidade
A variação espacial de uma variável regionalizada pode ser, dependendo do fenômeno, grande ou pequena, mas deve existir uma certa continuidade ponto a ponto; Anisotropias
A regionalização pode apresentar anisotropias quando apresenta variações graduais numa direção e rápida ou irregular em outra;
5.3 O variograma
O variograma é a ferramenta básica que permite descrever quantitativamente a variação no espaço de um fenômeno regionalizado (Huijbregts, 1975). A natureza estrutural de um conjunto de dados (assumido pela variável regionalizada) é definida a partir da comparação de valores tomados simultaneamente em dois pontos, segundo uma determinada direção. A função variograma 2γ(h) é definida como sendo a esperança matemática do quadrado da diferença entre os valores de pontos no espaço, separados por uma distância h, conforme a seguinte expressão: 2γ (h )
= E {[ Z ( x + h ) − Z ( x )]2 }
ou em termos computacionais: 2γ (h )
=
1
n
.∑ [ Z ( x + h) − Z ( x)]
2
n i =1
33
onde: 2γ(h) é a função variograma; n é o número de pares de pontos separados por uma distância h; Z(x) é o valor da variável regionalizada no ponto x; Z(x+h) é o valor da variável regionalizada no ponto (x+h). Comumente utiliza-se da função semivariograma, que é simplesmente a metade da função variograma: γ (h )
n
1
=
.∑ [ Z ( x + h) − Z ( x)]
2
2n i =1
Imaginava-se a expressão da função semivariograma como empírica, conforme proposta por Matheron (1971), mas Journel (1989) mostrou que ela não é nada mais que o momento de inércia medido num diagrama de dispersão entre os valores de Z(x+h) versus Z(x), como apresenta-se a seguir. Para o desenvolvimento da relação de dependência entre valores (x i,yi), separados por uma distância h, considere-se o diagrama de dispersão da Figura 15, apresentado por Journel (1989). A distância di entre o i-ésimo ponto (xi,yi) e a reta ideal é: d i
=
xi
− yi
cos 45 o
elevando a distância ao quadrado tem-se: 2
d i
=
1 2
( xi − yi )2
Figura 15: Representação do par de pontos (x i,yi) no diagrama de dispersão (Journel, 1989). Havendo n pares de pontos, pode-se calcular o momento de inércia em torno da reta de 45o, como: γ xy
=
1
n
.∑ d i2
n i =1
,
34
ou γ xy
=
1
n
.∑
1
n i =1 2
( xi − yi )2 =
1
n
.∑ ( xi
2n i=1
− y i )2
Quanto maior a dispersão, maior o momento de inércia e menor a correlação. Se não houver dispersão - todos os pares de pontos caem sobre a reta 45o - o momento de inércia é zero e o coeficiente de correlação é igual a 1 (máxima correlação). Como pode ser visto, tem-se uma medida eficiente da dependência espacial por meio do momento de inércia do conjunto de pontos separados por uma certa distância em relação à reta 45o. 5.4 Relação entr e semivariograma e a função covariância
Como na Estatística Clássica, pode-se definir a média e a variância de uma variável regionalizada, de acordo com as seguintes relações: m
= E [ Z ( x )]
Var [ Z ( x )]
= E [ Z ( x ) − m]2
A variância é conhecida em notação geoestatística como C(0), ou seja, a covariância para distância de separação nula. Da mesma forma, pode-se definir a covariância C(h), entre pontos separados por uma distância h: C (h )
= E [ Z ( x + h ). Z ( x )] − m 2
(4)
A função variograma 2γ(h) pode também ser expressa em termos de variância C(0) e da covariância C(h), de acordo com o seguinte desenvolvimento: 2
2γ (h ) = E [ Z ( x + h ) − Z ( x )] γ (h ) =
1 2
[
]
E Z ( x + h ) − 2 Z ( x + h ). Z ( x ) + Z ( x ) 2
2
aplicando-se as propriedades (b) e (c) da média (item 4.3.1) obtém-se: γ (h ) =
1
{ E [ Z ( x + h)]− 2 E [ Z ( x + h). Z ( x)] + E [ Z ( x)]} 2 2
2
(5)
desenvolvendo-se a expressão da variância tem-se: 35
2
C (0) = E [ Z ( x ) − m]
C (0) = E Z ( x ) − 2 Z ( x ).m + m 2
2
aplicando-se novamente as propriedades (b) e (c) da média (item 4.3.1): C (0) = E Z ( x ) 2
− 2mE [ Z ( x )] + m 2
como E [ Z ( x )] = m , tem-se: C (0) = E Z ( x ) 2
− 2mm + m 2 = E Z 2 ( x ) − m 2
ou E Z ( x ) 2
= C (0) + m 2
(6)
admitindo-se a estacionaridade, ou seja, que a média do quadrado da variável regionalizada no ponto (x) é igual àquela no ponto (x+h): E Z 2 ( x )
= E Z 2 ( x + h)
(7)
Substituindo-se (4), (6) e (7) em (5), a função γ(h) fica: γ (h ) =
1
γ (h) =
1
{C (0) + m 2 2
2
− 2[C (h) + m 2 ]+ C (0) + m 2 }
[2C (0) − 2C (h )]
portanto: γ (h ) = C (0) − C (h )
(8)
Como a função variograma é uma medida da variância das diferenças nos valores da variável regionalizada entre pontos separados por uma distância h, pontos mais próximos, por estarem correlacionados terão essa variância pequena, aumentando à medida que os pontos se distanciam. Ao contrário da função covariância, que é grande para distâncias pequenas diminuindo à medida que a distância aumenta, pois esta função mede a correlação entre pontos separados por uma distância h. A função variograma é usualmente representada sob a forma gráfica denominada variograma e a da função covariância é 36
denominada covariograma. A Figura 16 mostra a relação entre a função variograma e a função covariograma.
Figura 16: Relação entre as funções variograma e covariograma. O variograma é determinado segundo uma direção predefinida, portanto a função γ(h) é vetorial. Na prática faz-se variogramas segundo várias direções da jazida, justamente para se conhecer a estrutura da mineralização. Este procedimento é denominado "análise estrutural" na literatura (e.g. Huijbregts, 1975; Olea, 1994). 5.5 Propriedades do variograma
A interpretação do variograma permite obter parâmetros que descrevem o comportamento espacial das variáveis regionalizadas. As principais propriedades do variograma, que podem ser vistas na Figura 17, são:
Figura 17: Desenho mostrando um variograma típico e suas propriedades. Ampl itude
É a distância a partir da qual as amostras passam a ser independentes (Figura 17). Em outras palavras, a amplitude reflete o grau de homogeneização entre as amostras, ou seja, quanto maior for a amplitude maior será a homogeneidade 37
entre as amostras. Nesse sentido, conforme Matheron (1971), o variograma dá um significado preciso da noção tradicional de zona de influência. A amplitude (a) é a distância que separa o campo estruturado (amostras correlacionadas) do campo aleatório (amostras independentes); Patamar
É o valor de variância no qual o variograma estabiliza-se (no campo aleatório); Efeito pepita
É o valor da função variograma na origem (h=0). Teoricamente esse valor deveria ser zero, pois duas amostras tomadas no mesmo ponto (h=0) deveriam ter os mesmos valores; entretanto, quando não é assim, atribui-se, esta diferença, geralmente, a erros de amostragem e/ou análise. Como os erros analíticos são desprezíveis, com os equipamentos disponíveis atualmente, o efeito pepita é atribuído a erros de amostragem e/ou à variabilidade natural do depósito. O efeito pepita também é chamado de variância aleatória (Figura 17).; Variância espacial
É dada pela diferença entre a variância a priori e o efeito pepita (Figura 17); Zona de influência
Uma feição resultante da análise dos parâmetros do variograma experimental é a determinação da zona de influência, que é um fenômeno de transição caracterizado exclusivamente por modelos de variograma que possuem patamar e amplitude definidos. Portanto, qualquer valor de Z(x) estará correlacionado com outros valores Z(x+h) que estiverem dentro de um raio “a” de x. Esta correlação, ou a influência de um valor em outro, decresce conforme Z(x+h) aproxima-se de “a”. 5.6 Anisotropias
Os variogramas determinados ao longo de diferentes direções da jazida podem mostrar variações distintas, como exemplificados pela Figura 18. A anisotropia pode ser geométrica (Figura 18A), quando a amplitude varia conforme as direções, mas sob um patamar constante; zonal (Figura 18B) quando a amplitude permanece constante e o patamar varia de acordo com a direção; e, por fim, a anisotropia mista (Figura 18C) onde variam tanto a amplitude quanto o patamar, ou seja, quando as várias direções resultam em diferentes variogramas.
38
Figura 18: Anisotropias: geométrica (A), zonal (B) e mista (C). A anisotropia pode ser identificada facilmente através da confecção e análise de variogramas direcionais. Após ajuste dos modelos, anota-se os valores de amplitude e patamar, com os quais constrói-se uma rosácea, à qual ajusta-se uma elipse visando a definição precisa da direção de anisotropia, bem como a quantificação dos eixos de maior e menor elongação. Exemplo de anisotropia geométrica
Num depósito eólico a permeabilidade deve ter uma amplitude maior na direção do vento em relação à amplitude na direção perpendicular. Exemplo de aniso tropia zonal
O variograma de um furo de sonda vertical mostra uma patamar maior que na direção horizontal. 5.7 Comport amento próximo à origem
O grau de continuidade da mineralização é dado pelo comportamento do variograma próximo à origem. Assim, quanto a este comportamento podem ser descritos quatro tipos básicos, a saber: Parabólico
O variograma descreve uma curva parabólica próximo à origem (Figura 19A) e representa um alto grau de continuidade das amostras selecionadas. Este tipo pode ser exemplificado por um variograma construído a partir de dados de espessura de uma camada; Linear
Caracterizado por um comportamento linear na origem, ou seja, por uma tangente oblíqua à origem (Figura 19B), representando uma continuidade média das amostras. Entenda-se por continuidade média das amostras como sendo uma grande homogeneidade destas a pequenas distâncias e uma progressiva perda 39
de homogeneidade com o aumento da distância. Este comportamento é típico de muitos depósitos minerais metálicos; Efeito pepita
Este tipo apresenta uma descontinuidade na origem, (Figura 19C). Esta descontinuidade pode ser reflexo de dois fatores não mutuamente exclusivos erros de medida na amostragem e micro variabilidades; Efeito pepita puro
É um tipo extremo de comportamento do variograma próximo à origem (Figura 19D) e reflete a variação espacial de um fenômeno de transição, onde para um dado valor de patamar a amplitude terá um valor infinitesimalmente menor que as distâncias de observação (Journel & Huijbregts, 1978). O efeito pepita puro é um fenômeno de difícil ocorrência em mineralizações, porém ressalta-se que neste caso não se deve utilizar o método geoestatístico de interpolação. Burguess & Webster (1980) citam que o termo efeito pepita teve origem na mineração de ouro, onde a inclusão de uma pepita de ouro em uma pequena amostra de um testemunho de sondagem é um evento aleatório.
Figura 19: Graus de continuidade da mineralização expressos pelo comportamento do variograma na origem: alto grau de continuidade (A); média continuidade (B): efeito pepita (C) e efeito pepita puro (D), segundo Bubenicek & Haas, 1969).
5.8 Domínio do variograma
O domínio de definição do variograma é chamado campo geométrico, o qual implica que o variograma é válido dentro desse domínio e, portanto, se este for alterado, o variograma deve ser recalculado. 40
O campo geométrico, considerado geologicamente homogêneo, deveria ser até certo ponto intrínseco ou independente da posição nas características representando a variabilidade da variável regionalizada (Bubenicek & Haas, 1969). Quando esta hipótese é verificada, reconhece-se uma lei de dispersão única no campo mineralizado denominada "lei intrínseca", de acordo com aqueles autores. Entretanto, esta característica intrínseca não se mantém quando se move o campo para a zona de borda. Nesta parte a mineralização não existe e, por isso, o valor da variável regionalizada pode ser zero e o conceito de valor médio nesse campo de um dado tamanho tornar-se-ia insignificante, como está ilustrado na Figura 20. 5.9 Cálcul o de variogramas experimentais
A obtenção de variogramas representativos depende fundamentalmente do número de pares de pontos, para diferentes distâncias, encontrado numa determinada direção. Portanto, as direções devem ser especificadas para colherem o máximo de informações. Basicamente pode-se obter variogramas horizontais e verticais, sendo que para os primeiros deve-se ainda especificar direções e aberturas para pesquisa de pontos para fins de cálculo da função variograma. Os variogramas assim obtidos servem para identificar e determinar possíveis anisotropias.
Figura 20: Domínios de definição do variograma, segundo Bubenicek & Haas (1969). O campo geométrico coincide com o depósito e, neste caso, um variograma intrínseco pode ser obtido (A); o campo geométrico engloba parte do depósito e uma zona não mineralizada, fazendo com que o variograma seja dependente da posição e tamanho do campo, além de apresentar variabilidade maior que aquela verificada no variograma intrínseco (B); o campo geométrico é muito maior que o depósito e o variograma tende a zero quando o tamanho do campo aumenta. Contudo, é possível definir um variograma transitivo que é independente do campo que engloba o depósito (C).
No caso de variogramas horizontais para uma malha de amostragem quadrada, com direção da linha base E-W, especifica-se quatro direções iniciais de pesquisa: E-W, N45oE, N-S, e N45oW. Geralmente os variogramas horizontais são especificados segundo a orientação da linha base e daí a 45 o, 90o e 135o, no sentido anti-horário. 41
A consistência dos pontos do variograma experimental irá depender exclusivamente do número de pares de amostras. Para fins práticos Journel & Huijbregts (1978) recomendam utilizar, no mínimo, 30-50 pares de amostras para cada ponto do variograma experimental. Acrescenta-se ainda que se deve sempre observar o número de pares de pontos usado para o cálculo de γ(h) próximo à origem do variograma. Neste sentido, deve-se cuidar que no momento do ajuste do modelo de variograma os pontos que definirão o efeito pepita, se houver, deverão conter o maior número de pares de pontos possíveis e deve-se descartar aqueles com um número muito inferior àquele locado imediatamente após. A comparação entre valores de amostras separadas por uma distância h é direta se os pontos de dados estiverem distribuídos segundo uma malha regular. Entretanto, quando os pontos de dados estiverem dispersos, deve-se fazer a pesquisa de amostras situadas a uma distância h, dentro de uma janela de pesquisa. Esta janela é definida, ao longo da direção do variograma, por um ângulo e por uma distância de tolerância, conforme pode ser observado na Figura 21. O ângulo de tolerância pode ser limitado levando-se em consideração a distância “percorrida” ao longo da direção, ou seja, quando a tolerância angular é estabelecida forma-se um triângulo (2D) ou um cone (3D) em torno da direção preferencial. O problema é que se não houver uma limitação, a área do triângulo ou o volume do cone tendem a crescer indefinidamente, englobando maior número de pontos. Para evitar isso, define-se a largura máxima, isto é, estabelece-se uma distância a partir da qual o triângulo ou cone ficam limitados a essa faixa. N
i a c â n s s o r e l a T o o P d
Passo 4
s o a s P
a u r r g i m a a L á x M
Passo 3 Passo 2
Passo 1 Passo 0
Tolerância Angular
Direção
E
Figura 21: Desenho mostrando a direção do variograma, os passos, a tolerância angular, a largura máxima e a tolerância do passo (modificado de Pannatier, 1994). 5.10 Modelos teóri cos de variogramas
O variograma como ferramenta básica será utilizado para calcular os valores da função variograma, para uma dada distância, os quais são necessários para a organização do sistema de equações de krigagem. O variograma de 42
pontos, dito experimental, não serve para esse fim, porque há necessidade de interpolação e, invariavelmente, os pontos apresentar-se-ão com uma certa dispersão, principalmente para distâncias grandes, quando o número de pares de amostras vai diminuindo. Assim, surge a necessidade de ajustar uma função matemática que descreva continuamente a variabilidade ou correlação espacial existente nos dados. O ajuste de uma função matemática ao variograma experimental é denominado modelagem de variogramas. Esta modelagem é feita de maneira interativa, onde a partir dos parâmetros do variograma (modelo, efeito pepita, amplitude e patamar), o variograma teórico é desenhado juntamente com os pontos do variograma experimental, e se o ajuste não for satisfatório, novos parâmetros são fornecidos sucessivamente, até que o ajuste seja considerado satisfatório. Os modelos de variogramas mais comuns na natureza estão ilustrados na Figura 22, conforme as equações apresentadas a seguir. Exponencial Gauss Esférico
⎡ ⎛ ⎛ h ⎞ ⎞⎤ γ (h) = C o + C⎢1 − exp⎜⎜ − ⎜ ⎟ ⎟⎟⎥ ⎝ ⎝ a ⎠ ⎠⎦ ⎣ ⎡ ⎛ ⎛ h ⎞ 2 ⎞⎤ γ (h) = C o + C⎢1 − exp⎜ − ⎜ ⎟ ⎟⎥ ⎜ ⎝ a ⎠ ⎟⎥ ⎢⎣ ⎝ ⎠⎦ ⎡ 3 ⎛ h ⎞ 1 ⎛ h ⎞ 3 ⎤ γ (h) = C o + C⎢ ⎜ ⎟ − ⎜ ⎟ ⎥ para h < a ⎢⎣ 2 ⎝ a ⎠ 2 ⎝ a ⎠ ⎥⎦ γ (h) = C o + C para h ≥ a
Foram apresentados os três modelos mais comuns na natureza e que podem resolver a maioria dos problemas na modelagem da correlação espacial de fenômenos geológicos. Obviamente existem outros modelos, mas não se justifica introduzi-los num texto introdutório de geoestatística. Cabe ressaltar que em qualquer modelo teórico que apresente efeito pepita, bem como no efeito pepita puro, o valor de γ(h)=0 para h=0, pois o variograma é descontínuo na origem.
Figura 22: Modelos teóricos de variogramas mais comuns no estudo de fenômenos espaciais geológicos.
43
6 ESTIMATIVAS POR KRIGAGEM ORDINÁRIA
Após a análise geoestatística, na qual os variogramas experimentais foram calculados e os modelos teóricos foram ajustados, passa-se ao cálculo de estimativas pela técnica da krigagem ordinária. A krigagem ordinária tem como característica principal a precisão local das estimativas, mas com perda da precisão global devido ao efeito de suavização (suavização da variância e do variograma). Por outro lado, a técnica da simulação estocástica tem sido preferida para estudo da variabilidade, pois a variância de krigagem não proporciona uma medida precisa da incerteza associada à estimativa. A simulação estocástica reproduz tanto o histograma como o variograma, mas com significativa perda de precisão local. Infelizmente, de acordo com Olea (1999), as realizações estocásticas não estão livres de erros na representação da realidade e os erros, para qualquer realização, são maiores que aqueles da estimativa de krigagem. Esta é uma característica menos atrativa da simulação estocástica (Olea, 1999). Portanto, não há uma solução pronta para obtenção de uma única imagem representativa que compartilhe tanto a precisão global e local. Assim, ambas aproximações necessitam de uma correção para representar apropriadamente o fenômeno espacial através de uma única imagem representativa. Muitos autores têm preferido tentar corrigir o efeito de suavização da krigagem ordinária, ao invés das simulações estocásticas. Cabe lembrar, contudo, que a precisão local é muito mais importante que a precisão global, quando o problema for a estimativa de recursos naturais, a partir de dados de amostragem. Na realidade, todo o esforço na procura de alternativas para determinação da variabilidade foi justificado pela falta de uma medida da variância do erro. Assim, este autor (Yamamoto, 2000) propôs recentemente uma alternativa ao cálculo da variância do erro através da variância de interpolação que será visto adiante. Antes de passar a estimativa propriamente dita, a krigagem como qualquer outro método de interpolação requer a definição de certas condições de controle visando estimativas de qualidade. Tais condições são: definição da fronteira convexa e da vizinhança local. 6.1 Definição da fronteira convexa
As estimativas só podem ser feitas dentro do domínio dos pontos de dados, que pode ser aproximado através da sua fronteira convexa. A fronteira convexa pode ser definida como o polígono convexo de área mínima que engloba os pontos de dados. Cabe notar que a maioria dos programas de geoestatística não permite a definição da fronteira e, conseqüentemente, estimando pontos fora do domínio dos pontos amostrados, sem nenhum significado prático ou real. O uso de limites para interpolação de dados evita a interpretação de dados espúrios criados por extrapolação matemática (Yamamoto, 1997). Detalhes de algoritmos para determinação de fronteira podem ser encontrados em Yamamoto (1997). A título de ilustração a Figura 23 apresenta um conjunto de pontos de dados e a sua fronteira convexa. Observe-se que apenas os nós da malha regular que estão dentro da fronteira convexa podem ser estimados. 44
A)
B)
100
100
80
80
60
60
40
40
20
20
0 0
20
40
60
80
100
0 0
20
40
60
80
100
Figura 23: Conjunto de pontos de dados (A) e sua fronteira convexa com o desenho dos nós da malha regular pertencentes à mesma (B).
6.2 6.2 Defin Definição ição da vizinhança local
A krigagem ordinária faz uso da correlação espacial existente entre amostras, modelada pela função variograma. Isto significa que somente as amostras dentro de um raio de influência (igual à amplitude) poderão ser utilizadas para a estimativa do valor da variável de interesse em um ponto não amostrado. Assim, a krigagem ordinária é uma técnica essencialmente local, ao contrário da superfície de tendência que é global (todas as amostras são consideradas para o ajuste de uma superfície em toda a área de estudo). Sendo a krigagem uma técnica local de estimativa, deve-se estabelecer estratégias para localização e pesquisa das amostras vizinhas mais próximas do ponto a ser estimado. A localização e busca de n amostras de furos vizinhos para definição do subconjunto de amostras a ser utilizado na estimativa local é um passo importante da krigagem ordinária. Pois, dependendo do modo de pesquisa, diferentes subconjuntos de amostras poderão ser definidos e, portanto, resultados distintos poderão ser obtidos. A escolha das n amostras de furos vizinhos deve ser feita de tal modo que garanta uma boa amostragem espacial, o que implica em evitar subconjuntos com agrupamentos de pontos. Agrupamentos de pontos ocorrem preferencialmente em arranjos aleatórios e semi-regulares. Assim, torna-se necessário estabelecer critérios de seleção de amostras que garantam uma boa amostragem espacial e, conseqüentemente, evitem os agrupamentos de pontos. Existem basicamente três critérios que podem ser aplicados para a definição da vizinhança local: n pontos mais próximos, n/4 pontos mais próximos por quadrante e n/8 pontos mais próximos por octante. 45
Pontos mais mais próximo s
Considere-se, por exemplo, que o subconjunto de pontos seja definido pelos oito pontos mais próximos, em relação ao ponto a ser interpolado, para os arranjos aleatório e semi-regular, como ilustrado na Figura 24A e 24B, respectivamente. Na Figura 24A, pode-se observar que a pesquisa dos vizinhos próximos, sem nenhuma restrição quanto à localização dos mesmos, resulta no agrupamento de pontos no quadrante nordeste, em detrimento dos demais, enquanto o quadrante sudoeste nem sequer foi amostrado. No arranjo semiregular da Figura 24B, verifica-se que somente os pontos situados ao longo de uma linha de pesquisa serão amostrados, se nenhuma restrição for imposta, caracterizando também um agrupamento de pontos. Em nenhum caso a amostragem espacial foi representativa em termos da reprodução do gradiente dos dados. A)
B)
Figura 24: Localização dos oito pontos mais próximos para o arranjo aleatório (A), localização dos oito pontos mais próximos para o arranjo semi-regular (B), modificado de Harbaugh et al. (1977). Assim, para se evitar agrupamentos de pontos foram estabelecidos critérios de seleção de amostras baseados na subdivisão da região do ponto a ser estimado em quatro ou oito setores, denominados respectivamente quadrante e octante. Quadrante
Pelo critério dos quadrantes, a região do ponto a ser estimado é subdividida em quatro setores e os n/4 pontos mais próximos por quadrante são selecionados. Observe-se na Figura 25 que o critério dos quadrantes proporciona uma melhor amostragem espacial. Como se pode observar na Figura 6-3B, a aplicação do critério dos quadrantes provocou a amostragem de pontos em duas linhas adjacentes de pesquisa. 46
A)
B)
Figura 25: Seleção de duas amostras por quadrante, para o arranjo aleatório (A) e para o arranjo semi-regular (B), adaptado de Harbaugh et al. (1977).
Octante
Utilizando o critério dos octantes, a região do ponto a ser estimado é subdividida em oito setores, nos quais são escolhidos os n/8 pontos mais próximos por octante são selecionados. A Figura 26 apresenta os resultados da seleção pelo critério dos octantes mostrando uma melhor distribuição espacial dos pontos amostrados. A)
B)
Figura 26: Seleção de uma amostra por octante, para o arranjo aleatório (A) e para o arranjo semi-regular (B), adaptado de Harbaugh et al. (1977). Embora a seleção de amostras pelo critério dos octantes resulte numa melhor distribuição espacial, há, por outro lado, o inconveniente de amostras mais distantes serem selecionadas para a estimativa do ponto. Sem dúvida, deve 47
existir um compromisso entre a representatividade da amostragem e a distância máxima das amostras selecionadas. Quadrante sólido
Arranjos semi-regulares ocorrem também com informações de furos de sonda, onde a densidade de amostragem ao longo dos furos é sempre maior que entre os furos. A Figura 27 ilustra o caso da seleção de amostras de furos de sonda para interpolação de ponto ou bloco na jazida, sem impor nenhuma restrição. A Figura 28 mostra a mesma situação anterior, porém com restrição de localização por "quadrantes", ou seja, o equivalente para o caso tridimensional em o que cada quadrante representa um setor com ângulo sólido de 90 . Observe-se na Figura 27 que, no caso de avaliação de jazidas, o critério de seleção por setor é importante i mportante para evitar a superamostragem de um determinado furo, em relação aos demais. Na Figura 28, aplicando-se o critério de seleção por setor, pode-se verificar que os quatro furos de sonda foram amostrados, melhorando a representatividade da amostragem, principalmente quando se estiver avaliando um bloco de cubagem.
Figura 27: Localização de oito amostras de furos de sonda mais próximas ao centro do bloco.
Figura 28: Seleção de uma amostra de furo de sonda mais próxima por setor (octante tridimensional), em relação ao centro do bloco. 48
A pesquisa de amostras de furos vizinhos para interpolação de pontos ou centros de blocos para fins de avaliação de recursos, deve ser feita sempre com aplicação do critério de seleção em setores para garantir uma boa amostragem espacial. Caso contrário, há riscos de se estar avaliando blocos com base na média de amostras do furo mais próximo. Número de amostras de furos vizinhos
Escolhido o critério para a seleção de amostras de furos vizinhos, deve-se definir o número de amostras a ser utilizado para estimativa do valor de interesse em um ponto não amostrado. O número de amostras não deve ser excessivamente pequeno, com o risco da interpolação resultar em valor semelhante ou muito correlacionado ao do ponto mais próximo, e nem excessivamente grande, com o risco da interpolação resultar num valor bastante suavizado, perdendo a característica de interpolação local. Assim, pode-se definir 8 amostras, que se ajustam perfeitamente aos critérios de quadrante (2 amostras por quadrante) ou octante (1 amostra por octante) no plano, ou então ao critério de octante tridimensional. Entretanto, nem sempre a condição inicial de 8 amostras de furos vizinhos será satisfeita, principalmente na borda do corpo de minério. Nesses casos, deve-se relaxar a condição inicial para um mínimo de 3 ou 4 amostras, dependendo se a estimativa estiver sendo feita em 2D ou 3D, respectivamente.
6.3 Definição da malha regular
O último passo antes do cálculo de estimativas pela krigagem ordinária consiste na definição da malha regular em 2D ou 3D. Mas, por quê definir uma malha regular? Porque a malha regular proporciona áreas ou volumes de mesmo tamanho permitindo assim fazer uma comparação de resultados, bem como uma maior facilidade computacional para representação gráfica em mapas ou projeção em perspectiva. A malha regular pode ser definida em 2D ou 3D, dependendo dimensionalidade dos dados. A Figura 29 apresenta uma malha regular 2D, cujos nós pertencentes à fronteira convexa serão estimados pela técnica da krigagem ordinária. No caso da malha 2D, pode-se tanto estimar o nó, como uma área em torno do nó da malha regular. Daí a diferença entre krigagem pontual e de bloco como se verá adiante. Como os dados em 3D são geralmente ligados a resultados da pesquisa mineral, a malha regular em 3D é definida em termos de blocos de cubagem (não mais pontos). Os blocos de cubagem têm a forma geral de paralelepípedos e suas dimensões devem ser compatíveis com a densidade média de amostragem nas três direções. Ao conjunto de blocos de cubagem que compõem o depósito denomina-se modelo tridimensional de blocos (Figura 30). 49
100
80
60
40
20
0 0
20
40
60
80
100
Figura 29: Malha regular 2D com a representação dos nós pertencentes à fronteira convexa.
Figura 30: Modelo tridimensional de blocos de um depósito hipotético. A abertura ideal da malha regular, baseada na prática de avaliação de recursos, seria igual à metade do espaçamento médio entre os furos de sonda. Segundo Vallée & Côte (1992), a krigagem de blocos com dimensão muito menor que a metade da malha de amostragem deveria ser evitada, pois tais estimativas exibem extrema variabilidade. Cabe ressaltar que no caso da malha regular 3D, a fronteira convexa é definida também para todos os níveis do modelo tridimensional de blocos.
50
6.4 Krigagem ordinária
Segundo Brooker (1979), as técnicas geoestatísticas de estimativa, baseadas no estudo da variabilidade espacial do corpo de minério, são superiores porque permitem o cálculo do erro associado às estimativas, chamado variância de krigagem. Ainda conforme o mesmo autor, a krigagem é o procedimento que permite calcular os ponderadores para uma dada configuração (bloco X disposição das amostras no espaço), com mínima variância de krigagem. A krigagem é feita após a conclusão dos estudos geoestatísticos, os quais poderão inclusive indicar a não aplicação deste método se o comportamento da variável regionalizada for totalmente aleatório. Os estudos geoestatísticos levam a definição de um modelo de variograma, que servirá para inferir os valores da função variograma ou covariograma que serão utilizados pelos métodos geoestatísticos de interpolação. Equações de krig agem
A krigagem é um método que permite estimar o valor desconhecido Z ∗ (x o ) associado a um ponto, área ou volume, a partir de um conjunto de n dados {Z(x i), i=1,n} disponíveis. O estimador Z ∗ (x o ) poderá ser obtido como uma combinação linear dos dados disponíveis, conforme: n
Z (x o ) = ∑ λ i .Z(x i ) ∗
(9)
i=1
Os ponderadores (λi, i=1, n) são obtidos da resolução de um sistema linear de equações, denominado sistema de equações de krigagem, conforme o desenvolvimento matemático. Para que o estimador Z ∗ (x o ) não seja enviesado, segundo Journel & Huijbregts (1978), basta garantir que: E Z(x o ) − Z ∗ (x o ) = 0 fazendo E[Z(x o )] = m e tendo que:
⎡n ⎤ n E[Z ∗ (x o )] = E⎢∑ λ i Z(x i )⎥ = ∑ λ iE[Z(x i )] ⎣ i=1 ⎦ i=1 n
E[Z (x o )] = m∑ λ i ∗
i=1
51
assim, a condição de não enviesamento para Z ∗ (x o ) fica: n
∑λ
i
=1
(10)
i=1
Como toda técnica de estimativa, a krigagem procura fazê-la com mínima variância. A variância do erro da krigagem é dada pela equação a seguir:
σ E2 = Var {Z(x o ) − Z ∗ (x o )} Expandindo a variância do erro, de acordo com Isaaks & Srivastava (1989), tem-se:
σ E2 = Cov{Z(x o )Z(x o )} − 2Cov Z ∗ (x o )Z(x o )}+ Cov Z ∗ (x o )Z ∗ (x o )}
(11)
Desenvolvendo cada termo do lado direito de (11), conforme Isaaks & Srivastava (1989), tem-se: Cov{Z(x o )Z(x o )} = Var {Z(x o )} = C(0 )
⎧⎡ ⎫ ⎤ 2Cov{Z ∗ (x o )Z(x o )} = 2Cov ⎨⎢∑ λ i Z(x i )⎥ Z(x o )⎬ ⎦ ⎩⎣ i ⎭ ⎧ ⎫ ⎧ ⎫ = 2E⎨∑ λ i Z(x i )Z(x o )⎬ − 2E⎨∑ λ i Z(x i )⎬E{Z(x o )} ⎩ i ⎭ ⎩ i ⎭ = 2∑ λ iE{Z(x i )Z(x o )} − 2∑ λ iE{Z(x i )}E{Z (x o )} i
i
= 2∑ λ i [E{Z(x i )Z(x o )} − E{Z(x i )}E{Z (x o )}] i
= 2∑ λ iC(x o − x i ) i
Cov{Z ∗ (x o )Z ∗ (x o )} = Var {Z ∗ (x o )} ⎧ ⎫ = Var ⎨∑ λ i Z(x i )⎬ ⎩ i ⎭ = ∑∑ λ i λ j C(x i − x j ) i
j
Assim, a expressão (11) torna-se:
σ E2 = C(0 ) − 2∑ λ iC(x o − x i ) + ∑∑ λ i λ jC(x i − x j ) i
i
j
52
O objetivo da krigagem é buscar o melhor conjunto de ponderadores, de tal modo que a variância do erro seja a mínima possível. Trata-se, portanto, de encontrar o mínimo da função variância do erro. Entretanto, como tal função tem n variáveis, o ponto de mínimo poderá ser determinado após aplicação da técnica dos multiplicadores de Lagrange (Converse, 1970), conforme colocação do problema a seguir: - minimizar a função:
σ E2 = C(0 ) − 2∑ λ i C(x o − x i ) + ∑∑ λ i λ j C(x i − x j ) i
i
j
- restrito a:
∑λ
j
j
= 1 ou
∑λ
j
−1= 0
j
Forma-se o lagrangiano:
⎛ ⎞ L(λ 1, λ 2 ,K, λ n , µ ) = C(0 ) − 2∑ λ i C(x o − x i ) + ∑∑ λ i λ j C(x i − x j ) − 2µ⎜⎜ ∑ λ j − 1⎟⎟ i i j ⎝ j ⎠ onde: L(λ 1, λ 2 ,K, λ n , µ ) é o lagrangiano; µ é o multiplicador de Lagrange. Para minimizar o lagrangiano, faz-se cada uma das derivadas parciais dL/dλi iguais a zero: dL = −2C(x o − x i ) + 2∑ λ j C(x i − x j ) − 2µ = 0 dλ i j
para i=1,n
e fazendo dL/d µ igual a zero: dL = ∑λ −1= 0 dµ j j Assim, a minimização da variância do erro, sujeita à condição de não enviesamento, resulta nas de equações de krigagem ou sistema de krigagem:
⎧∑ λ j C(x i − x j ) − µ = C(x o − x i ) para i = 1,n ⎪ j ⎨ ⎪∑ λ j = 1 ⎩ j
(12)
53
Em termos matriciais, as equações de krigagem são representadas como segue:
⎡ C ( x1 − x1 ) ⎢C ( x − x ) ⎢ 2 1 ⎢ M ⎢ ⎢C ( xn − x1 ) ⎢⎣ 1
C ( x1 − x2 ) L
− x2 ) L
C ( x2
M
L
− x2 ) L
C ( xn
1
L
C ( x1 − x n )
⎡ λ 1 ⎤ ⎡ C ( xo − x1 )⎤ ⎥ ⎢ λ ⎥ ⎢C ( x − x )⎥ C ( x2 − xn ) 1 ⎥ ⎢ 2⎥ ⎢ o 2 ⎥ ⎥ M M⎥ ⋅ ⎢ M ⎥ = ⎢ M ⎥ ⎢ ⎥ ⎢ ⎥ C ( x n − x n ) 1⎥ ⎢ λ n ⎥ ⎢C ( xo − x n )⎥ ⎥⎦ 1 0⎥⎦ ⎢⎣− µ ⎥⎦ ⎢⎣ 1 1⎤
O sistema de equações de krigagem também pode ser escrito em termos da função semivariograma, visto que esta e a função covariância estão relacionadas conforme a equação (8):
⎧∑ λ j γ (x i − x j ) + µ = γ (x o − x i ) para i = 1,n ⎪ j ⎨ ⎪∑ λ j = 1 ⎩ j
(13)
Em notação matricial, o sistema de equações de krigagem escrito em termos da função semivariograma torna-se:
⎡γ ( x1 − x1 ) ⎢γ ( x − x ) ⎢ 2 1 ⎢ M ⎢ ⎢γ ( xn − x1 ) ⎢⎣ 1
γ ( x1 − x2 ) L γ ( x1 − x n ) γ ( x 2
− x2 ) L M
γ ( xn
L
− x2 ) L 1
L
⎡ λ 1 ⎤ ⎡γ ( xo − x1 )⎤ ⎥ ⎢λ ⎥ ⎢γ ( x − x )⎥ γ ( x2 − xn ) 1 ⎥ ⎢ 2⎥ ⎢ o 2 ⎥ ⎥ M M⎥ ⋅ ⎢ M ⎥ = ⎢ M ⎥ ⎢ ⎥ ⎢ ⎥ γ ( x n − xn ) 1⎥ ⎢λ n ⎥ ⎢γ ( xo − x n )⎥ ⎥⎦ 1 0⎥⎦ ⎢⎣ µ ⎥⎦ ⎢⎣ 1 1⎤
Variância de krigagem
A minimização da variância do erro resulta na variância de estimativa ou de krigagem ordinária, conforme segue: 2 σ KO = C(0 ) − 2∑ λ i C(x o − x i ) + ∑∑ λ i λ j C(x i − x j ) i
o termo
∑∑ λ λ C(x i j
i
i
i
j
− x j ) pode ser derivado do primeiro conjunto de equações
j
(12), conforme Isaaks & Srivastava (1989):
∑ λ C(x j
i
− x j ) − µ = C(x o − x i )
para i=1,n
j
escrevendo as n equações para i=1,n e somando: 54
∑ ∑ λ λ C(x i j
i
− x j ) − µ∑ λ i = ∑ λ i C(x o − x i )
i
j
e lembrando que
i
∑λ
i
i
= 1, tem-se:
i
∑ ∑ λ λ C(x i j
i
− x j ) = ∑ λ i C(x o − x i ) + µ
i
j
i
substituindo este resultado na expressão da variância de krigagem ordinária, temse: 2 σ KO = C(0 ) − ∑ λ i C(x o − x i ) + µ
(14)
i
A variância de krigagem em termos da função semivariograma torna-se: 2 σ KO = ∑ λ i γ (x o − x i ) + µ
(15)
i
Variância de interpolação
A variância de krigagem é homoscedástica, ou seja, ela é independente dos valores dos pontos de dados usando para obter o estimador Z ∗ ( xo ) (Olea, 1991). A Figura 31 ilustra porque a variância de krigagem não é uma medida completa da incerteza.
Figura 31: Estimativa de blocos a partir da mesma configuração de dados (segundo Armstrong, 1994). Observe-se que em (A) a variância do erro deveria ser menor que em (B) por causa de dados mais consistentes. Como os arranjos de dados são idênticos em ambos os casos, as variâncias de krigagem são idênticas e também as estimativas de krigagem se não há anisotropia (Armstrong, 1994). A variância de interpolação fará uso tanto 55
do conjunto de pesos { λi, i=1,n} e dos valores dos dados z(x i), medindo a maior variabilidade existente no bloco B. Journel & Rossi (1989) já haviam concluído que a variância de krigagem mede apenas a configuração espacial dos dados; como a variância de krigagem depende apenas do variograma que é global, ela independente dos valores locais dos pontos de dados. Assim, Yamamoto (2000) propôs uma expressão para o cálculo da variância de krigagem, que é determinada como a média ponderada das diferenças ao quadrado entre os valores dos pontos de dados e a estimativa Z ∗ ( xo ) , como segue: 2 S o
n
= ∑ λ i [ Z ( xi ) − Z ∗ ( xo )]
2
(16)
i =1
Segundo Yamamoto (2000), a variância de interpolação apresenta as seguintes propriedades: o
o o
Corresponde à propriedade de exatidão da krigagem ordinária, isto é, se o ponto a ser estimado coincide com um ponto de dado, então o peso deste ponto é igual a um com todos os outros pesos iguais a zero, portanto S o2 = 0 ; É proporcional à dispersão dos pontos de dados; Usa indiretamente a distância estrutural do variograma através do peso da krigagem ordinária λ i . Quanto mais influente o ponto de dado maior o seu peso.
Essa expressão foi introduzida por Yamamoto (1989) para definir uma variância de interpolação associada a teores estimados através das equações multiquádricas em depósitos minerais. Yamamoto (1991) estendeu a definição para calcular a variância de interpolação usando os pesos da krigagem ordinária. Correção de pesos negativo s
A expressão (16) pode resultar em variância negativa quando alguns pesos da krigagem ordinária forem negativos, devido à existência de agrupamentos de pontos de dados associada a uma distribuição estatística com forte assimetria positiva (distribuições lognormais). Da mesma forma que teores negativos devem ser evitados, pesos negativos também devem ser evitados. Assim, surgiram na literatura várias propostas para eliminação de pesos negativos da krigagem ordinária. Entre as propostas existentes, o autor tem adotado aquela proposta por Journel & Rao (1996), que consiste em adicionar o módulo do maior peso negativo (c) a todos os pesos: τ i
=
λ i + c
∑n j=1 (λ j + c ) 56
Os pesos corrigidos podem ser agora substituídos nas equações (9) e (16) para cálculo do valor estimado e da variância de interpolação, respectivamente. Distribuição de probabilidade
Como os pesos da krigagem ordinária são todos positivos e têm uma soma igual a um, eles podem ser interpretados como probabilidades condicionais associadas aos n pontos de dados locais (Journel & Rao, 1996). Em qualquer localização dada xo, classifique os n pontos de dados vizinhos em ordem crescente: z(x 1 ) ≤ z(x 2 ) ≤ L ≤ z(x n )
Os pesos da krigagem ordinária associados aos valores z(x 1), z(x2), ... , z(xn) da variável aleatória Z constituem a distribuição de probabilidade de Z. Os pares ordenados (z(xi), λi) constituem a função de probabilidade de Z. A função de distribuição acumulada condicional é então modelada como: α
F(x o , z α ) = ∑ λ i , i =1
ou seja, pode-se determinar a probabilidade da variável aleatória z ser menor ou igual a zα no ponto xo, condicional aos n pontos de dados mais próximos de x o. Tipos de krigagem
As equações de krigagem permitem determinar o conjunto de ponderadores {λi, i=1,n} associados ao conjunto de dados disponíveis {Z(xi),i=1,n}, que combinados conforme a equação (9), resulta na estimativa do valor desconhecido Z ∗ (x o ) . Conforme o domínio que se estima, tem-se: - krigagem pontual; - krigagem de bloco. Krigagem pontual
A krigagem pontual tem por objetivo estimar uma localização não amostrada. Vamos considerar situação apresentada na Figura 32, para ilustrar a krigagem pontual.
57
1
2
4 3
Figura 32: Configuração dos pontos de amostragem (+) para krigagem pontual na localização do ponto marcado com círculo. Os dados dos pontos de amostragem, selecionados pelo critério dos quadrantes, encontram-se na Tabela 8. Tabela 8: Coordenadas dos pontos de amostragem. Ponto Coord. X Coord. Y Variável 0* 250 200 ¿ 1 300 350 56 2 125 250 46 3 200 125 37 4 400 150 42 * Ponto a ser estimado pela krigagem ordinária. Seja considerado válido para os dados em estudo, o seguinte modelo de variograma: 3 ⎡ h h ⎞ ⎤ ⎛ γ (h ) = 5 + 15⎢1.5 − 0.5⎜ ⎟ ⎥ ⎝ a ⎠ ⎥⎦ ⎢⎣ a γ (h ) = 20
para h
< 800
para h ≥ 800
Para estimativa no ponto 0, precisamos dos ponderadores da krigagem ordinária que aplicados na expressão (9) resultará no valor estimado. Assim, para determinar os ponderadores precisamos resolver o sistema de equações de krigagem. Nessa demonstração vamos utilizar o sistema de equações escrito em termos da função semivariograma (13). A organização do sistema de equações (13) começa com o cálculo da matriz dos termos γ xi − x j , conforme segue:
Para o cálculo da função semivariograma, por exemplo, entre as amostras 1 e 2, determina-se inicialmente a distância entre as mesmas: 58
(300 − 125)2 + (350 − 250)2 = 201,56
d ( x1 , x 2 ) =
A distância encontrada é convertida em função semivariograma: 3 ⎡ 201.56 201.56 ⎞ ⎤ ⎛ γ ( x1 − x 2 ) = 5 + 15⎢1.5 − 0.5⎜ ⎟ ⎥ = 10.55 ; 800 800 ⎝ ⎠ ⎢⎣ ⎥⎦
Para os elementos da diagonal principal, observe-se que as distâncias são iguais a zero e, portanto, γ ( xi − xi ) = 0 ;
Repete-se o procedimento para todos os pares de amostras, obtendo-se a matriz dos termos γ xi − x j :
⎡ 0 ⎢10.55 ⎢ ⎢11.71 ⎢ ⎢11.13 ⎢⎣ 1
10.55 11.71 11.13 1⎤ 0
9.05
12.86
9.05
0
10.55
12.86 10.55 1
1
Resolvendo-se o sistema obtêm-se os ponderadores da krigagem ordinária:
[λ 1 = 0.21
1
0
⎡ λ 1 ⎤ ⎡9.39⎤ ⎥ ⎢λ ⎥ ⎢8.75⎥ 1 ⎥ ⎢ 2⎥ ⎢ ⎥ 1⎥ ∗ ⎢λ 3 ⎥ = ⎢7.52⎥ ⎥ ⎢ ⎥ ⎢ ⎥ 1⎥ ⎢λ 4 ⎥ ⎢9.39⎥ 0⎥⎦ ⎢⎣ µ ⎥⎦ ⎢⎣ 1 ⎥⎦
λ 2
= 0.231
λ 3
= 0.33
λ 4
= 0.23
= 0.58]
O valor da variável no ponto 0 é estimado como: Z ∗ ( xo ) = 0.21 * 56 + 0.23 * 46 + 0.33 * 37 + 0.23 * 42 = 44.21
O valor da variável de interesse no ponto 0 (não amostrado) é igual a 44.21, cujo resultado é bastante razoável em relação às amostras fornecidas. Krigagem de bloco
A krigagem de bloco é uma técnica de estimativa do teor médio em painéis ou blocos de cubagem, tratando-se, portanto de uma técnica desenvolvida exclusivamente para mineração. A estimativa de painéis ou blocos é muito diferente da estimativa pontual, à medida que áreas ou volumes maiores devem ser representados pelos pontos de amostragem. Assim, é certo que apenas a estimativa de um único ponto no centro daquelas unidades não será suficiente para representá-las, devendo haver uma diferença composicional entre o ponto estimado e a unidade lavrada. À essa 59
diferença composicional denomina-se erro de estimativa que depende fundamentalmente da amostragem, e esta por sua vez é função da variabilidade natural do depósito que está se avaliando. Certamente o erro de estimativa associado a krigagem de bloco será menor que aquele associado a krigagem pontual (ponto calculado no centro do bloco para representá-lo). O princípio da krigagem de bloco é baseado na subdivisão do bloco de cubagem em sub-blocos, os quais são avaliados individualmente e compostos para o bloco original, conforme o Teorema da Combinação das Estimativas de Krigagem (Journel & Huijbregts, 1978). Este teorema prova que tanto as estimativas como os ponderadores dos sub-blocos individuais podem ser combinados para dar origem à estimativa ou ponderadores médios do bloco de cubagem. Da mesma forma os vetores dos valores da função semivariograma, entre as amostras e os centros dos sub-blocos, podem ser combinados para dar origem ao vetor médio, dos valores da função semivariograma entre amostras e o bloco. A subdivisão do bloco de cubagem deve ser feita dentro dos limites máximos recomendados por Journel & Huijbregts (1978), conforme se reproduz na Tabela 9. Tabela 9: Limites máximos recomendados para a subdivisão do bloco a ser estimado (Journel & Huijbregts, 1978). dimensão do domínio 1 2 3
número de pontos 10 6x6 4x4x4
Estabelecendo um bloco de 100 x 100 e adotando uma subdivisão mínima do bloco (2 x 2), tem-se a configuração apresentada na Figura 33. 1
2
4 3
Figura 33: Configuração do bloco e centros de sub-blocos para avaliação pela krigagem ordinária de bloco. As coordenadas dos centros dos sub-blocos encontram-se na Tabela 10. 60
Para o cálculo dos ponderadores da krigagem ordinária vamos resolver o mesmo sistema de equações de krigagem (13) utilizado para a krigagem pontual, porém com a diferença do vetor do lado direito do sistema que levará em consideração a subdivisão em sub-blocos. Na verdade, a matriz dos termos γ xi − x j é a mesma para a krigagem pontual, pois as amostras são exatamente as mesmas. Assim, precisamos calcular o vetor da função semivariograma entre as amostras e os sub-blocos como mostra a Figura 34. Tabela 10: Coordenadas dos centros dos sub-blocos. Subbloco sb1 sb2 sb3 sb4
1
A
Leste (m) 275 225 225 275
Norte (m) 225 225 175 175
1
B
2
2
4
3
1
C
1
D
2
4
3
2
3
4
3
4
Figura 34: Esquema mostrando o cálculo da função semivariograma entre a amostra 1 e todos os sub-blocos (A); para a amostra 2 (B); para a amostra 3 (C) e para a amostra 4 (D). 61
A seguir tem-se as etapas para cálculo do vetor médio da função semivariograma entre as amostras e os sub-blocos:
Para ilustrar o procedimento, considere-se a amostra 1 (Figura 34A) em relação aos sub-blocos {sbi, i=1,4}, cujas distâncias são assim calculadas: d ( x sb1 − x1 ) = d ( x sb 2
− x1 ) = (300 − 225)2 + (350 − 225)2 = 145.77
d ( x sb3 − x1 ) = d ( x sb 4
(300 − 275)2 + (350 − 225)2 = 79.06
(300 − 225)2 + (350 − 175)2 = 190.39
− x1 ) = (300 − 275)2 + (350 − 175)2 = 176.78
As distâncias encontradas são convertidas em valores da função semivariograma: 3 ⎡ 79.06 79.06 ⎞ ⎤ ⎛ γ ( x sb1 − x1 ) = 5 + 15⎢1.5 − 0.5⎜ ⎟ ⎥ = 7.22 800 800 ⎝ ⎠ ⎥⎦ ⎢⎣ 3 ⎡ 145.77 145.77 ⎞ ⎤ ⎛ γ ( x sb 2 − x1 ) = 5 + 15⎢1.5 − 0.5⎜ ⎟ ⎥ = 9.05 800 800 ⎝ ⎠ ⎣⎢ ⎦⎥ 3 ⎡ 190.39 190.39 ⎞ ⎤ ⎛ γ ( x sb3 − x1 ) = 5 + 15⎢1.5 − 0.5⎜ ⎟ ⎥ = 10.25 800 800 ⎝ ⎠ ⎥⎦ ⎢⎣ 3 ⎡ 176.78 176.78 ⎞ ⎤ ⎛ − 0.5⎜ γ ( x sb 4 − x1 ) = 5 + 15⎢1.5 ⎟ ⎥ = 9.89 800 800 ⎝ ⎠ ⎢⎣ ⎥⎦
O valor médio da função semivariograma entre a amostra 1 e o bloco, pode ser calculado:
γ ( xo
− x1 ) =
7.22 + 9.05 + 10.25 + 9.89 4
= 9.1025
Dessa forma são calculados os demais valores médios:
62
γ ( xo
− x 2 ) =
9.23 + 7.88 + 8.49 + 9.65 4
= 8.8125
γ ( xo
− x 3 ) =
8.49 + 7.88 + 6.57 + 7.52 4
= 7.6150
γ ( xo
− x 4 ) =
9.05 + 10.25 + 9.89 + 8.56 4
O sistema de equações de krigagem de bloco torna-se:
⎡ 0 ⎢10.55 ⎢ ⎢11.71 ⎢ ⎢11.13 ⎢⎣ 1
10.55 11.71 11.13 1⎤ 0
9.05
12.86
9.05
0
10.55
12.86 10.55 1
1
0 1
⎡ λ 1 ⎤ ⎡9.1025 ⎤ ⎥ ⎢λ ⎥ ⎢8.8125 ⎥ 1 ⎥ ⎢ 2⎥ ⎢ ⎥ 1⎥ ∗ ⎢λ 3 ⎥ = ⎢7.6150 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ 1⎥ ⎢λ 4 ⎥ ⎢9.4375 ⎥ 0⎥⎦ ⎢⎣ µ ⎥⎦ ⎢⎣ 1 ⎥⎦
Resolvendo-se o sistema obtém-se os ponderadores da krigagem de bloco:
[λ 1 = 0.24
= 9.4375
λ 2
= 0.22
λ 3
= 0.32
λ 4
= 0.22
µ = 0.55]
O valor da variável no ponto 0 é estimado como: Z ∗ ( xo ) = 0.24 * 56 + 0.22 * 46 + 0.32 * 37 + 0.22 * 42 = 44.64
O valor médio do bloco foi ligeiramente superior ao valor do ponto, devido à pequena variabilidade dos dados, bem como à dimensão do bloco. A fim de mostrar na prática os limites máximos de subdivisão de blocos, reproduz-se na Tabela 11, as estimativas obtidas para vários níveis de subdivisão, segundo Isaaks & Srivastava (1989). Tabela 11: Exemplos de estimativa do valor médio de V (parâmetro de interesse) dentro de blocos de 10 x 10 m 2 usando a krigagem ordinária de bloco e várias malhas de discretização dentro do bloco (segundo Isaaks & Srivastava, 1989). Centro do bloco E N 80 80 100 80 80 90 100 90 80 100 100 100 80 110 100 110
1x1 584,67 408,53 538,36 460,13 497,66 530,37 781,17 591,13
malhas de subdivisão 2x2 4x4 6x6 576,41 574,30 573,98 418,29 419,19 419,38 519,89 520,58 520,53 479,73 480,35 480,52 547,87 549,40 550,13 513,32 513,56 513,47 737,04 732,29 731,06 580,73 578,75 578,74
10 x 10 573,81 419,47 520,47 480,61 550,51 513,42 730,41 578,72 63
6.5 Validação cruzada
A validação cruzada consiste em estimar a localização de um ponto de dado eliminando-se o valor do mesmo do conjunto de pontos de dados. Repetindo-se este procedimento um número de vezes igual ao número de pontos do conjunto, tem-se ao final o valor estimado e o valor verdadeiro, além das variâncias de krigagem e de interpolação. Trata-se da técnica do jackknife herdada da estatística clássica para predição do erro de estimativa. Uma alternativa à validação cruzada seria a simulação, porém tem a grande desvantagem de produzir realizações baseadas em populações hipotéticas e não na população amostrada como faz a validação cruzada. A Figura 35A apresenta um conjunto de pontos de dados que será submetido à validação cruzada. Na Figura 35B tem-se a indicação de um ponto de dado sendo estimado a partir dos 8 vizinhos mais próximos escolhidos pelo critério dos quadrantes. Como a estimativa é feita pela krigagem ordinária pontual, os resultados da validação cruzada podem ser usados para aferir a modelagem de semivariograma. O resultado típico da validação cruzada é apresentado na forma de uma diagrama de dispersão dos valores da validação cruzada em função dos valores reais, como ilustra a Figura 36. O resultado ideal de uma validação cruzada seria a reta de regressão estar mais próxima da bissetriz e que a dispersão em torno desta reta fosse mínima. Observe-se que a validação cruzada pode ser utilizada para aferir a modelagem do semivariograma. Pode-se assim testar aproximações diferentes, mas se após procedimentos de tentativa e erro, os pontos da validação cruzada não apresentarem um bom ajuste em torno da bissetriz, significa que há um enviesamento condicional, causado pelo efeito de suavização das estimativas por krigagem ordinária.
A
B
Figura 35: Conjunto de pontos de dados (A) e validação cruzada de um ponto do conjunto, o qual será estimado a partir dos 8 pontos vizinhos mais próximos selecionados pelo critério dos quadrantes (B). 64
15
A Número de dados = 224 D Coeficiente de c orrelação = 0.725 A Correlação rankeada = 0.698 Z U R C O Ã Ç A D 10 I L A V
5
0 0
5
10
15 VA LOR REAL
Figura 36: Resultado típico de uma validação cruzada Pode-se observar na Figura 36 que os maiores valores estão em geral subestimados, enquanto os menores valores estão superestimados. Há, portanto, nos resultados da validação cruzada o efeito de suavização da krigagem ordinária. A reta de regressão confirma essa observação. Efeito d e suavização
Estimativas baseadas na fórmula da média ponderada, tal como a krigagem ordinária (9), apresentam uma variabilidade reduzida, que é referida na literatura como efeito de suavização, de acordo com Isaaks & Srivastava (1989, p.420). Além disso, conforme esses autores, a utilização de um maior número de amostras tende geralmente a aumentar a suavidade das estimativas. Como uma conseqüência do efeito de suavização, pequenos valores são geralmente superestimados enquanto valores altos são subestimados, caracterizando um enviesamento condicional do estimador resultante (Goovaerts, 1997, p. 370). O enviesamento condicional dado que o valor estimado z ∗ ( x ) é maior que um valor de referência zc é, segundo Journel et al. (2000): E Z ∗ ( x ) − Z ( x ) Z ∗ ( x ) > z c
≠0
O efeito de suavização surge como um sério problema na detecção de padrões de valores extremos do atributo, tais como zonas de alta permeabilidade ou zonas ricas em metal (Goovaerts, 1997, p.370). Além disso, o efeito de suavização é desigual no espaço, sendo zero nos pontos de dados e aumentando à medida que a localização x distancia-se dos pontos de dados (Journel et al., 2000, p. 791).
65
O efeito de suavização é devido a um déficit na variância do estimador da krigagem e para a krigagem simples (KS) é precisamente igual à variância de krigagem como mostrado por Journel et al. (2000, p. 791): ∗ ( x ) Var { Z ( x )}− Var Z KS
2 = σ KS ( x )
Similarmente, o déficit de variância para a krigagem ordinária (KO) pode ser derivado da seguinte expressão (Yamamoto, 2000, p. 507): ∗ ( x ) Var { Z ( x )}− Var Z KO
2 = σ KO (x ) + 2µ ≥ 0
onde µ é o multiplicador de Lagrange. Yamamoto (2000, p. 493) também mostrou que o déficit de variância poder ser expresso em termos da variância de interpolação como: ∗ ( x ) Var { Z ( x )}− Var Z KO
= E S o2 ≥ 0
que pode ser interpretado como a suavização da variância de interpolação do ∗ ( x ) . Esta expressão faz sentido como uma estimador da krigagem ordinária Z KO medida do déficit de variância porque quão maior é a variância de interpolação S o2 , tomada em média sobre todos os valores de dados possíveis, maior a suavização da variância de interpolação do estimador da krigagem (Yamamoto, 2000, p. 493). Como conseqüência do efeito de suavização, as estimativas pela krigagem ordinária não reproduzem o histograma e o variograma. A Figura 37 mostra comparativamente o histograma original dos pontos de dados (Fig. 37A) e o histograma dos valores estimados (Fig. 37B). O efeito de suavização pode ser observado nos coeficientes de variação que passa de 0,746 dos valores originais para 0,513 para os valores estimados. A distribuição dos valores estimados apresenta, conseqüentemente, menor assimetria em relação à distribuição original. Na Figura 38 tem-se os variogramas da imagem de referência (a partir da qual foram amostrados os pontos de dados originais) e da imagem estimada pela krigagem ordinária. É notável o patamar menor da imagem krigada, dado pelo efeito de suavização (menor variância).
66
A)
B) 15
25
%
N úm ero de dados = Média = Desvio padrão = Coeficiente de variação = Máximo = Quar til super ior = Mediana = Quartil inferior = Mínimo =
10
%
224 3.137 2.339 0.746 14.577 4.191 2.735 1.469 0.290
N úm er o de dados = 2015 Média = 3.147 Desvio padrão = 1.615 Coeficiente de variação = 0.513 Máximo = 14.577 Quartil superior = 3.828 Mediana = 2.914 Quartil inferior = 2.092 Mínimo = 0.290
20 15 10
5
5
0
0
5
10
15 VARIÁVEL
0 0
5
10
15
ESTIMATIVA (KO)
Figura 37: Histograma da distribuição dos dados originais (A) e dos dados estimados pela krigagem ordinária (B). 10
M A R G O 8 I R A V I M E 6 S
4
2
0 0
200
400
600
800
1000 DISTANCE
Figura 38: Semivariogramas para a imagem de referência (símbolos cheios) e para a imagem estimada pela krigagem ordinária (símbolos vazios). Legenda: quadrado = N135o e círculo = N45o. 6.6 Classific ação de recurso s/reservas minerais
A variância de estimativa ou variância de krigagem foi proposta como uma medida da incerteza associada à estimativa feita por meio da krigagem ordinária. Contudo, como demonstrado por diversos autores (Journel, 1986, Olea, 1991, entre outros), a variância de krigagem mede apenas a configuração espacial dos pontos de dados e, por isso, não reconhece a dispersão local dos mesmos. A dispersão local é de importância vital para fins de classificação de reservas, pois discrimina regiões de alta e baixa variabilidade, as quais podem ter atribuídas classes de baixa e alta confiabilidade, respectivamente. A variância de krigagem, ao utilizar a dispersão média do depósito medida por meio do variograma, não consegue discriminar regiões de alta e baixa variabilidade. Por outro lado, a
67
variância de interpolação é uma alternativa que mede a dispersão local, bem como reconhece o efeito proporcional quando existente (distribuições lognormais). Diversos autores (Froidevaux, 1982, Isaacks & Srivastava, 1989, Vallée & Côte, 1992 e Wober & Morgan, 1993) têm proposto calcular o erro com base no intervalo clássico em torno da estimativa, usando a distribuição normal e o desvio padrão da krigagem. Nestas aproximações, a distribuição normal não é adequada, pois o tamanho das amostras em problemas de avaliação de reservas é sempre inferior a 40-60 amostras/bloco avaliado. Diehl & David (1982) e Wellmer (1983) propuseram, também com base no desvio padrão da krigagem, utilizar o Teorema do Limite Central para calcular o erro em torno de uma estimativa. De fato, o Teorema do Limite Central proporciona uma boa aproximação para uma distribuição normal. Contudo, estes autores não consideraram o intervalo de confiança para determinação do erro, mas simplesmente o coeficiente de variação multiplicado pelo valor crítico de t, como segue: ERRO =
σ E Z ∗ ( xo )
t gl ,ns
onde σE é o desvio padrão de estimativa no sentido genérico podendo ser substituído pelo desvio padrão de krigagem ou de interpolação, e t gl,ns é o valor crítico da distribuição t de Student para gl graus de liberdade e nível de significância ns. Yamamoto & Rocha (1996) propuseram a seguinte expressão para o cálculo do erro de estimativa ou tolerância permitida, para fins de classificação de reservas minerais: ERRO =
1 IC 2 NC (%) Z ∗ ( xo )
=
σ E .t gl ,ns ∗
Z ( xo ). n
100(%)
(17)
onde 1 2 IC NC (%) é a metade do intervalo de confiança a um nível de confiança (NC), σE é o desvio padrão da estimativa, t gl,ns é o valor crítico de t para gl graus de liberdade e um nível de significância ns, Z ∗ ( xo ) é o valor estimado e n o número de amostras utilizadas para a estimativa. A Figura 39 mostra a curva da distribuição t centrada na estimativa X e o seu intervalo de confiança a um determinado nível de confiança. A distribuição t só depende do número de graus de liberdade. Quando o número de graus de liberdade tende ao ∞, a distribuição t passa para uma normal. Contudo, a partir de 40-60 amostras, a distribuição t passa praticamente para uma normal. Como dificilmente utiliza-se 40-60 amostras para avaliação de um bloco de cubagem, o erro de uma estimativa deve ser determinado segundo o Teorema do Limite Central (expressão 17) e, portanto, considerando a distribuição.
68
Figura 39: Curva da distribuição t de Student, utilizada para o cálculo do intervalo de confiança da estimativa X a um determinado nível de confiança NC. A expressão (17) tem sentido para cálculo do erro associado a estimativas pontuais, pois tem n em seu denominador, sendo n o número de amostras. Contudo, em estimativas de bloco, onde o bloco é subdividido em nsb sub-blocos, a n deveria ser substituída pela nsb , pois, na verdade, são utilizados nsb subblocos para a estimativa de teor e da variância de interpolação. Assim, a expressão para cálculo do erro associado à estimativa de bloco fica: ERRO =
σ E .t gl ,ns Z ∗ ( xo ). nsb
100(%)
(18)
O modelo para classificação de recursos/reservas minerais, proposto para adoção pelo DNPM (1992), foi baseado no modelo australiano atualizado em 1996 pela AIMM- Australasian Institute of Mining and Metallurgy (AIMM, 1996). Neste modelo, os recursos minerais são classificados, de acordo com nível crescente de conhecimento do depósito, em: inferido, indicado e medido. A partir da consideração de fatores econômicos, de mineração, metalúrgicos, mercadológicos, legais, ambientais, sociais e políticos, os recursos indicado e medido poderão vir a transformar-se em reservas provável e provada, respectivamente (AIMM, 1996). Os recursos minerais que apresentam indicação de economicidade, segundo DNPM (1992), deveriam ser classificados na fase de pesquisa mineral, para suportar os estudos de viabilidade técnico-econômica que seguem. A classificação de recursos minerais, conforme a proposta do DNPM (1992), foi baseada no modelo da ONU ( apud Valente, 1980). Tais modelos classificam os recursos minerais, com uma confiabilidade de 95%, em medido, quando os erros são menores que 20%, indicado quando os erros estão entre 20 e 50% e inferido quando os erros são superiores a 50%. Yamamoto & Rocha (1996) consideraram que o nível de confiança igual a 95%, para fins de classificação de recursos minerais não é adequado, pois seria aplicável somente em depósitos minerais extremamente homogêneos. Um nível de confiança igual a 90% foi considerado razoável por Yamamoto & Rocha (1996), contudo, consideraram também que poderiam ser aplicados somente a alguns depósitos com baixa a média 69
variabilidade. Segundo Koch & Link (1971), em aplicações envolvendo avaliação e outros assuntos de economia mineral, o uso de um nível de confiança de 90% é consistente com os altos riscos associados à indústria mineral, riscos estes que se tornam claros com as altas taxas de retorno da mineração em relação às indústrias de transformação. Assim, Yamamoto & Conde (1999) propuseram a adoção dos níveis de erro propostos pela ONU (apud Valente, 1980) e DNPM (1992), sob um nível de confiança igual a 90%, para fins de classificação de recursos minerais. Este modelo, em comparação com aqueles propostos por Diehl & David (1982), Wellmer (1983) e DNPM (1992), encontra-se na Tabela 12. Os erros adotados no modelo proposto são calculados pela expressão (23), com a substituição do desvio padrão de estimativa pelo desvio padrão de interpolação. Pelos motivos expostos anteriormente, o desvio padrão de interpolação apresenta maior precisão e confiabilidade para fins de classificação de reservas, como se procurará demonstrar neste artigo. Como se sabe, a variabilidade dos depósitos minerais é fortemente afetada por uma componente denominada variabilidade natural. A variabilidade natural resulta da interação dos processos geológicos, os quais afetam a distribuição espacial dos teores e a morfologia do depósito. Assim, dependendo da variabilidade natural do depósito em estudo, as classes de erros propostas na Tabela 12 podem não estar adequadas e, portanto, não permitirão classificar os recursos em medido, indicado e inferido, conforme esperado. Muitas vezes, dependendo do nível de detalhe da pesquisa mineral, os recursos ainda não podem ser classificados como medidos, devido à tolerância máxima de 20%, e sim como indicados. Nestas situações, a certeza geológica, expressa pela confiabilidade dos dados e pela continuidade da mineralização, deveria prevalecer sobre os erros de tolerância. Tabela 12: Modelo para classificação de recursos minerais, em comparação com as propostas de Diehl & David (1982), Wellmer (1983) e DNPM (1992). Recursos⇒ Reservas⇒ Diehl & David (1982) Wellmer (1983) ONU e DNPM (1992) Yamamoto & Conde (1999)
Medido Provada Provável Erro: ±10% Erro: ±20% N.C.: > 80% N.C.:60-80% Erro: ±10% Erro: ±20% N.C.: 90% N.C.: 90% Erro: ±0-20% N.C.: 95% Erro: 0-20% N.C.: 90%
Indicado Possível Erro: ±40% N.C.:40-60% Erro: ±30% N.C.: 90% Erro: ±20-50% N.C.: 95% Erro: 20-50% N.C.: 90%
Inferido Inferida Erro: ±60% N.C.:20-40% Erro: ±50% N.C.: 90% Erro: >±50% N.C.: 95% Erro: > 50% N.C.: 90%
N.C. = é o nível de confiança conforme a Figura 4.
70
7 ESTIMATIVAS POR COKRIGAGEM ORDINÁRIA
7.1 Definições Básicas de Isotop ia e Heterotopi a
Segundo Wackernagel (1998), medidas de diferentes variáveis em um dado domínio podem localizar-se tanto em um mesmo ponto de amostragem como em pontos diferentes. Assim, dependendo da coincidência ser total ou parcial, define-se os termos isotopia e heterotopia. A heterotopia pode ser parcial quando em alguns pontos de amostragem foram medidas as variáveis em estudo, enquanto nos demais pontos de amostragem uma ou outra variável. A Figura 40 ilustra em (A) o caso de isotopia, onde as duas variáveis de interesse foram analisadas em todos os pontos de amostragem, em (B) a heterotopia é parcial, pois como pode ser observado, em apenas alguns pontos de amostragem as duas variáveis de interesse foram analisadas e, por fim, em (C) exemplifica-se a heterotopia total onde não existem pontos de amostragem comuns às duas variáveis. (A) Ag
Au 0 8
Ag
Au
Ag
Ag
0 4
Ag
0
Ag
20
Au
40
60
Ag
0
Au
0 2
Au
Ag
80
Au
Au
0 2
Au
Ag
0 6
Ag
Ag
Au
0 2
0
Ag Au
0 4
Ag
Au
Ag
0 6
Ag Au
0 8
Ag Au
Au
Ag
Au
0 8
Ag
0 6
0 4
(C )
Ag
Au
Au Au
Au
(B)
Ag
Au
0
20
Au
40
60
80
Ag
Ag 0
0
20
Au
40
60
80
Figura 40: (A) exemplifica um caso onde a amostragem é isotópica, ou seja, as duas variáveis de interesse encontram-se analisadas em todos os pontos de amostragem. (B) apresenta um caso de heterotopia parcial onde as variáveis de interesse possuem apenas alguns pontos de amostragem onde a analise foi realizada em ambas e (C) representa o caso de heterotopia total onde não há pontos de amostragem comuns para ambas variáveis. 7.2 O variograma cruzado
O variograma cruzado foi definido por Matheron em 1965 como a generalização natural do variograma e pode ser escrito como: γ ij (h ) =
1 2
[
]
E {[ Z i ( x + h ) − Z i ( x )]* Z j ( x + h ) − Z j ( x ) }
71
para funções aleatórias intrínsecas multivariadas (de ordem zero) deve satisfazer (Chilès & Delfiner 1999):
⎧ E [ Z i ( x + h ) − Z i ( x )] = 0 ⎨Cov[ Z ( x + h ) − Z ( x ), Z ( x + h ) − Z ( x )] = 2γ (h ) i i j j ij ⎩
para i
= 1,L, p
existe e depende apenas de h
Ainda segundo Chilès & Delfiner (1999), o variograma cruzado possui duas vantagens sobre o covariograma cruzado, a saber: 1. não assume variâncias finitas e; 2. a estimativa do variograma cruzado não é contaminada pela estimativa das médias. Quando existe, a relação entre o variograma cruzado e a covariância cruzada é: γ ij (h ) = C ij (0) −
1 2
[C (h) + C (− h)] + 1 [C (h) − C (− h)] ij
ij
2
ij
ij
Wackernagel (1998) descreve a o variograma cruzado como uma função obviamente ímpar que satisfaz a desigualdade: γ ii (h )γ jj (h ) ≥ γ ij (h )
2
pois o quadrado da covariância dos incrementos de duas variáveis é limitado ao produto das variâncias do incremento correspondente. Em outras palavras, a covariância obtida no variograma cruzado não pode ser maior que o produto das variâncias dos variogramas diretos. Assim como o variograma experimental, o variograma cruzado, é uma função direcional e diferencia-se daquele por considerar duas variáveis diferentes no cálculo da variância. Os outros procedimentos seguem exatamente aqueles para o cálculo de variogramas experimentais. Observando-se a equação do variograma cruzado pode-se ver que é a expansão da equação do variograma: 2γ (h ) = E [ Z ( x + h ) − Z ( x )]
2
= E [ Z ( x + h ) − Z ( x ) * Z ( x + h ) − Z ( x)] ,
da qual diferencia-se por ser o segundo termo do segundo membro uma variável diferente daquela analisada no primeiro termo. Algumas questões práticas devem ser mencionadas para o variograma cruzado. a. uma dada estrutura (por exemplo: efeito pepita) não pode estar presente no variograma cruzado ( γ 12 ) se não estiver presente nos variogramas diretos ( γ 1 ou γ 2 ) ;
72
b. uma dada estrutura (por exemplo: efeito pepita) pode estar presente nos variogramas diretos ( γ 1 ou γ 2 ) e não estar presente no variograma cruzado ( γ 12 ) . 7.3 O Modelo Li near de Corregionalização
Segundo Wackernagel (1998), a regionalização multivariada de um conjunto de funções aleatórias pode ser representado por um modelo linear espacial multivariado. O modelo linear de corregionalização é a soma de modelos de covariâncias proporcionais. Em notação matricial C (h ) = C ij (h ) é a matriz covariância com dimensões pxp e similarmente Γ(h ) = γ ij (h ) é a matriz variância, esse modelo assume a forma (Chilès & Delfiner 1999): C (h ) =
n
∑ B C (h) k
1
k
n
ou Γ(h ) = ∑ Bk γ k (h ) 1
Do modelo linear de corregionalização resulta que as variâncias espaciais das estruturas que compõem os modelos de variograma cruzados são dependentes das variâncias espaciais dos modelos de variograma diretos. Destaca-se que, para o modelo linear de corregionalização ser honrado, o quadrado da variância espacial do variograma cruzado deve ser menor ou igual que o produto das variâncias espaciais dos variogramas diretos, por exemplo: Tomem-se os seguintes modelos de variograma ajustados aos variogramas experimentais: h ⎞ ⎛ h ⎞⎤ e γ = 5 + ⎡9sph⎛ h ⎞⎤ onde γ é o variograma da = 3 + ⎡⎢5sph⎛ ⎜ ⎟ + 12sph⎜ ⎟ ⎜ ⎟ 2 1 ⎢⎣ ⎝ 50 ⎠ ⎝ 135 ⎠⎥⎦ ⎝ 135 ⎠⎥⎦ ⎣ variável primária e γ 2 é o variograma da variável secundaria. O modelo teórico do variograma cruzado será admissível se possuir efeito pepita máximo e variância espacial máxima iguais a:
γ 1
h ⎞⎤ = 3,873 + ⎡⎢10,392 sph⎛ ⎜ ⎟⎥ observa-se ainda que o variograma cruzado não 135 ⎝ ⎠⎦ ⎣ apresenta a estrutura com amplitude igual a 50m, pois o variograma direto da variável secundária não a apresenta.
γ 12
7.4 Cokri gagem ordinária
A cokrigagem é a extensão natural da krigagem quando dados multivariados e variogramas multivariados ou modelos de covariância podem ser calculados (Wackernagel, 1998). O procedimento da cokrigagem consiste em 73
estimar uma variável de interesse em um ponto específico com base nas informações vizinhas da própria variável e nas informações disponíveis para variáveis auxiliares ou secundárias. Segundo Olea (1999), a cokrigagem pode ser descrita como um procedimento de estimativa verdadeiramente multivariado, pois o modelo lida com dois ou mais atributos (variáveis) em um mesmo domínio. Em geoestatística quando duas ou mais variáveis regionalizadas são definidas em um campo aleatório são chamadas de corregionalização. No caso da cokrigagem, quando existe ausência de análise de uma variável em um determinado ponto de amostragem, esta ausência não interfere, ou enviesa, os resultados. Pelo contrário, a cokrigagem apresenta sua melhor performance quando se verifica esta situação. Assim como na krigagem, existem vários algoritmos de cokrigagem, dentre os quais se destacam a cokrigagem ordinária, a cokrigagem simples e cokrigagem co-localizada. Tratar-se-á, aqui, apenas da cokrigagem ordinária, uma vez que este é o algoritmo mais utilizado por não requerer que a média populacional seja conhecida como no caso da cokrigagem simples. A cokrigagem co-localizada necessita que em todos os pontos de amostragem, ou na maioria deles, a variável secundária seja conhecida, com informação apenas parcial da variável primária. Equações de cok rigagem ordinária
A estimativa da cokrigagem ordinária é feita através de uma combinação linear de pesos λ pi , a partir de dados de diferentes variáveis localizados em pontos de amostragem na vizinhança de um ponto x0 (Wackernagel, 1998). Ainda segundo este autor, cada variável é definida em um conjunto de amostras de tamanho n p e o estimador é definido como: N
* p0
Z
n p
( x0 ) = ∑∑ λ pi Z p ( xi ) p =1 i =1
onde o índice p0 refere-se a uma variável específica de um conjunto de N variáveis. O número de amostras n p depende do índice p das variáveis. Segundo Wackernagel (1998), no arcabouço da hipótese intrínseca conjunta deseja-se estimar uma variável específica em um conjunto de N variáveis com base em um erro de estimativa, o qual deve ser nulo em média. Esta condição, segundo o mesmo autor, é satisfeita determinando-se pesos cuja soma seja um para a variável de interesse (primária) e seja zero para a variável auxiliar (secundária), conforme: n p
∑ λ
p i
i =1
⎧1 se p = p0 = δ pp = ⎨ ⎩0 se p ≠ p 0 0
74
Expandindo-se a expressão do erro médio de estimativa, conforme desenvolvimento de Wackernagel (1998), tem-se: E [ Z
* p0
( x0 ) − Z p
⎡ n ⎤ n N N n ⎢ ⎥ ( x0 )] = E ⎢∑∑ λ pi Z p ( xi ) − ∑ λ pi Z p ( x0 ) − ∑∑ λ pi Z p ( x0 )⎥ = p = 0 i =1 i =1 ⎢ p =1 i =1 ⎥ 1 23 123 p ≠ p ⎣ ⎦ 1 0 p
p
p
0
0
0
0
N
n p
= ∑∑ λ pi E [ Z p ( xi ) − Z p ( x0 )] = 0 144 4 244 4 3
p =1 i =1
0
Continuando, a variância do erro de estimativa fica: 2 N n ⎡ ⎤ ⎛ ⎞ σ E 2 = E ⎢⎜⎜ ∑∑ λ pi Z p ( xi ) − Z p ( x0 )⎟⎟ ⎥ ⎢⎣⎝ p =1 i =1 ⎠ ⎥⎦ p
0
⎧− 1 se p = p0 Introduzindo-se os pesos λ p0 = −δ pp = ⎨ que estão incluídos nos 0 ≠ se p p ⎩ 0 somatórios, pode-se reduzir a expressão da variância de estimativa para: 0
2 N n ⎡ ⎤ ⎛ ⎞ p 2 ⎜ ⎟ ⎢ σ E = E ⎜ ∑∑ λ i Z p ( xi ) ⎟ ⎥ ⎢⎣⎝ p =1 i =0 ⎠ ⎥⎦ p
Inserindo-se variáveis aleatórias fictícias posicionadas na origem, pode-se formar incrementos:
Z p (0 )
arbitrariamente
2 ⎡⎛ ⎛ ⎤ ⎞ ⎞ n n ⎜ ⎟ N ⎜ ⎟ ⎢ ⎥ σ E 2 = E ⎢⎜ ∑ ⎜ ∑ λ pi Z p ( xi ) − Z p (0 )∑ λ pi ⎟ ⎟ ⎥ i =1 ⎢⎜⎜ p =1 ⎜⎜ i =0 1 23 ⎟⎟ ⎟⎟ ⎥ ⎢⎣⎝ ⎝ ⎠ ⎠ ⎥⎦ 0 p
p
2 ⎡⎛ N n ⎞ ⎤ = E ⎢⎜⎜ ∑∑ λ pi ( Z p ( xi ) − Z p (0 ))⎟⎟ ⎥ ⎢ p =1 i =0 144 244 3 ⎥ ⎠ ⎥⎦ incrementos ⎢⎣⎝ p
P ( xi , x j ), que não é Definido-se a covariância cruzada dos incrementos C pq invariante à translação, tem-se:
N
2 E
σ
N
n p
nq
= ∑∑∑∑ λ pi λ q j C pqP ( xi , x j ) p =1 q =1 i = 0 j = 0
75
Para converter as covariâncias dos incrementos para variogramas, deve-se assumir que as covariâncias cruzadas dos incrementos são simétricas. Com esta hipótese obtém-se o valor da translação invariante, como: N
2 E
σ
n p
N
N
n p
nq
= 2∑∑ λ γ pp ( xi − x0 ) − γ p p ( x0 − x0 ) − ∑∑∑∑ λ pi λ q j γ pq ( xi − x j ) p i
p =1 i =1
0
0
0
p =1 q =1 i =1 j =1
Após a minimização, na qual as restrições dos pesos geraram N multiplicadores de Lagrange p , ter-se-á o sistema de cokrigagem ordinária:
⎧ N n q ⎪∑∑ λ j γ pq ( xi − x j ) + µ p = γ pp ( xi − x0 ) para p =1,L, N ; i = 1,L, n p ⎪ q =1 j =1 ⎨n ⎪ ∑ λ p j = δ pp para p = 1, L, N ⎪⎩ j =1 q
0
p
0
em termos matriciais será escrito, em sua forma reduzida e para uma variável primária e uma secundária, como:
⎡ ⎢ C ⎢ 11 ⎢ ⎢ ⎢ ⎢ C 21 ⎢ ⎢1 1 ⎢ ⎢⎣0 0
1 1
C 12
1 0 C 22
0
0
1
0
0
0 0
0
1 1 1
0
⎡ λ 11 ⎤ ⎡ ⎤ ⎢ 1⎥ 0⎥ ⎢ λ 2 ⎥ ⎢C 01 ⎥ ⎥ 1 ⎢ ⎥ ⎥ 0⎥ ⎢ λ 3 ⎥ ⎢ ⎢ ⎥ 2 ⎢ ⎥ 1⎥ λ 1 ⎥ ⋅ ⎢ 2 ⎥ = ⎢ ⎥ onde 1⎥ ⎢ λ 2 ⎥ ⎢C 02 ⎥ ⎢ ⎥ ⎢ ⎥ 1⎥ λ 32 ⎥ ⎢ ⎥ ⎢ ⎥ 0⎥ ⎢ µ 1 ⎥ ⎢ 1 ⎥ 0⎥⎦ ⎢⎣ µ 2 ⎥⎦ ⎢⎣ 0 ⎥⎦ 0⎤
C pq
é uma matriz covariância 3x3
e a variância de cokrigagem será: N
2 CKO
σ
n p
= ∑∑ λ pi γ pp ( xi − x0 ) + µ p − γ p p ( x0 − x0 ) p =1 i =1
0
0
0
0
Observa-se que o sistema de equações de cokrigagem pode ser escrito em termos de variogramas e, para tal, deve-se apenas inverter o sinal dos multiplicadores de Lagrange. REFERÊNCIAS BIBLIOGRÁFICAS
Aitchison, J. & Brown, J. A. C. 1957. The lognormal distribution, with special reference to its uses in economics. Apud Koch & Link (1971). 76
Armstrong, M., 1994, Is research in mining geostats as dead as a dodo ? in: Dimitrakopoulos, R., ed., Geostatistics for the next century: Kluwer Academic Publishers, Dordrecht, p. 303-312. Australasian Institute of Mining and Metallurgy-AIMM. 1996. Australasian code for reporting of identified mineral resources and ore reserves. The AusIMM Bulletin (4):1-19. Barnes, M.P. 1980. Computer assited mineral appraisal and feasibility. New York, Soc. Min. Eng. 167p. Blais, R. A. & Carlier, P. A. 1968. Applications of geostatistics in ore valuation. apud Olea, R.A. 1975. Optimum mapping techniques using regionalized variable theory. Kansas, Kansas Geol. Survey. 137p. (Series on Spatial analysis, 2). Bon, E. H. 1979. Exploration techniques employed in The Pulau Tujuh tin discovery. Trans. Inst. Min. Metall., 88(A13-A22). Brooker, P.I. 1979. Kriging. Eng. & Min. J., 180(9):148-153. Bubeniceck, L. & Haas, A. 1969. Method of calculation of the iron ore reserves in the Lorraine deposit. In: The decade of digital computing in the mineral industry. The American Institute of Mining, Metallurgical and Petroleum Engineers, Inc. New York. p.179-210. Burguess, T. M. & Webster, R. 1980. Optimal interpolation and isarithmic mapping of soil properties, I: the semivariogram and punctual kriging. J. Soil Sci., 31:333-341. Chilès, J-P & Delfiner, P. 1999, Geostatistics: Modeling Spatial Uncertainty: John Wiley & Sons, Inc. New York, 1 st. Ed., 695p. Cochran, W.G. 1963. Técnicas de amostragem. Rio de Janeiro, Fundo de Cultura. 555p. Converse, A.O. 1970. Otimização. Trad. de A. Rios Neto. São Paulo, EDART/ EDUSP. 262p. Davis, J. C. 1986. Statistics and data analysis in geology. New York, John Wiley. 550p. Departamento Nacional da Produção Mineral-DNPM, 1992. Bases técnicas de um sistema de quantificação do patrimônio mineral brasileiro. Brasília, DNPM. 28p. (Relatório do Grupo de Trabalho, Portaria n o 03 de 24/10/90, da Secretaria Nacional de Minas e Metalurgia). 77
Diehl, P. & David, M. 1982. Classification of ore reserves/resources based on geostatistical methods. CIM Bull., 75(838):127-136. Exell, R.H.B. 1998. Practical mathematics, random variables. World Wide Web, http://www.crosswinds.net/~exell//. Fonseca, J. S. & Martins, G. A. 1982. Curso de estatística. São Paulo, Editora Atlas. 286p. Froidevaux, R. 1982. Geostatistics and ore reserve classification. CIM Bull. 75(843):77-83. Goovaerts, P., 1997, Geostatistics for natural resources evaluation: Oxford University Press, New York, 483p. Gy, P. 1968. Theory and practice of sampling broken ores. Montreal, Can. Inst. Min. Metall., Special Vol. 9, p.5-10. Harbaugh, J.W.; Doveton, J.H.; Davis, J.C. 1977. Probability methods in oil exploration. New York, John Wiley. 269p. Huijbregts, C. J. 1975. Regionalized variables and quantitative analysis of spatial data. In: Davis, J.C. & McCullagh, M.J. (ed) Display and anaysis of spatial data. New York, John Wiley. p.38-53. Instituto de Pesquisas Tecnológicas do Est. de São Paulo - IPT. 1989. Estudo e implementação de novas técnicas e desenvolvimento de programas em geoestatística: aplicações em geologia. São Paulo, IPT. 98p. (IPT Relatório 27.799). Isaaks, E.H. & Srivastava, R.M. 1989. An introduction to applied geostatistics. New York, Oxford University Press. 561p. Journel, A.G. 1986. Geostatistics: models and tools for the Earth Sciences. Math. Geol., 18(1):119-140. Journel, A. G. 1989. Fundamentals of geostatistics in five lessons. Washington, American Geophysical Union. 40p. Journel, A. G. & Huijbregts, C.J. 1978. Mining geostatistics. London, Academic Press. 600p. Journel, A.G., and Rao, S.E., 1996, Deriving conditional distributions from ordinary kriging: Stanford Center for Reservoir Forecasting, Stanford, 25p. (Report #9) Journel, A.G., and Rossi, M.E., 1989, When do we need a trend model in kriging ?: Math. Geology, v. 21, no. 7, p. 715-739. 78
Journel, A., Kyriakidis, P.C., and Mao, S., 2000, Correcting the smoothing effect of estimators: a spectral postprocessor: Math. Geology, v. 32, no. 7, p. 787-813. Kim, Y. C. 1990. Introductory geostatistics and mine planning. Tucson, The University of Arizona. 139p. Koch Jr., G.S. & Link, R.F. 1971. Statistical analysis of geological data. New York, Dover Publishing. Vol. I. 375p. Krige, D. G. 1978. Lognormal-De Wijsian geostatistics for ore evaluation. Johannesburg, South Afr. Inst. Min. Metall. 50p. Landim, P. M. B. - 1985. Introdução a Geoestatística. Laboratório de Geomatemática. Publicação Didática n.º 3. Instituto de Geociências e Ciências Exatas, Campus de Rio Claro. UNESP. 69p. Lapin, L. L. 1982. Statistics for modern business decisions. Apud Davis, J. C. (1986). Matheron, G. 1965. Les variables régionalisées et leur estimation. Apud Huijbregts, C.J. 1975. Regionalized variables and quantitative analysis of spatial data. Matheron, G. 1971. The theory of regionalized variables and its applications. Fontainebleau, École Nationale Supérieux das Mines de Paris. 211p. (Les Cahiers de Centre de Morphologie Mathématique de Fontainebleau, 5). Miesch, A.T. 1967. Theory of errors in geochemical data. Washington, USGS. (USGS Prof. Paper, 574-A). 17p. Olea, R.A., 1991, Geostatistical glossary and multilingual dictionary: Oxford University Press, New York, 175p. Olea, R. 1994. Fundamentals of semivariograms estimation, modeling and usage. In: Yarus, J.M. & Chambers, R.L. (eds) Stochastic modeling and geostatistics – principles, methods and case studies. Tulsa, AAPG (AAPG Computer Application in Geology no 3). P.27-36. Olea, R., 1999, Geostatistics for engineers and earth scientists: Kluwer Academic Publishers, Boston, 303p. Pannatier, Y. - 1994. MS-WINDOWS Programs for Exploratory Variography and Variogram Modeling in 2D. In “Statistics of Spatial Processes: Theory and Applications. Bary. Italy. Royle, A. G. 1979. Why geostatistics ? Eng. & Min. J., 180(5):92-101. 79
Spiegel, M. R. - 1967. Estatística. São Paulo, McGraw Hill. 580p. Valente, J. 1980. Avaliação geoestatística com parametrização da jazida fosfática de Anitápolis. In: CONGR. BRAS. GEOL., 31, Camboriú, 1980. Anais...Camboriú, SBG, v.2, p.1243-1257. Vallée, M. & Côte, D. 1992. The guide to the evaluation of gold deposits: integrating deposit evaluation and reserve inventory practices. CIM Bull., 85(957):50-61. Wackernagel, H. 1998. Multivariate Geostatistics: Springer-Verlag, Berlin, 2 nd. Ed.,291p. Waeny, J.C. de C. 1979. Recomendações para nomenclatura de confiabilidade metrológica. São Paulo, IPT. (Publicação Interna, 2). 23p. Wellmer, F.W. 1983. Classification of ore reserves by geoestatistical methods. ERZMETALL, 36(7/8):315-321. Wober, H.H., and Morgan, P.J., 1993, Classification of ore reserves based on geostatistical and economic parameters. CIM Bulletin, v. 86, no. 966, p. 7376. Yamamoto, J.K. 1997. Convex_Hull – A Pascal program for determining the convex hull for planar sets. Computers & Geosciences, v. 23, no. 7, p. 725738. Yamamoto, J.K., 2000, An alternative measure of the reliability of ordinary kriging estimates: Math. Geology, v. 32, no. 4, p. 489-509. Yamamoto, J.K. & Rocha, M.M. 1996. Revisão e recomendações para o cálculo e classificação de reservas minerais. Rev. Bras. Geoc., 26(4):243-254. Yamamoto, J.K. & Conde, R.P. 1999, Classificação de reservas minerais usando a variância de interpolação: Rev. Bras. Geoc., 29(3):465-474.
80
ANEXO 1: DISTRIBUIÇÃO NORMAL PARA X ENTRE 0 E 3,49 E AS INTEGRAIS Q(X) CORRESPONDENTES. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4
0.00 5000 4602 4207 3821 3446 3085 2743 2420 2119 1841 1587 1357 1151 0968 0808 0668 0548 0446 0359 0287 0228 0179 0139 0107 0082 0062 0047 0035 0026 0019 0013 0010 0007 0005 0003
0.01 4960 4562 4168 3783 3409 3050 2709 2389 2090 1814 1562 1335 1131 0951 0793 0655 0537 0436 0351 0281 0222 0174 0136 0104 0080 0060 0045 0034 0025 0018 0013 0009 0007 0005 0003
0.02 4920 4522 4129 3745 3372 3015 2676 2358 2061 1788 1539 1314 1112 0934 0778 0643 0526 0427 0344 0274 0217 0170 0132 0102 0078 0059 0044 0033 0024 0018 0013 0009 0006 0005 0003
0.03 4880 4483 4090 3707 3336 2981 2643 2327 2033 1762 1515 1292 1093 0918 0764 0630 0516 0418 0336 0268 0212 0166 0129 0099 0075 0057 0043 0032 0023 0017 0012 0009 0006 0004 0003
0.04 4840 4443 4052 3669 3300 2946 2611 2296 2005 1736 1492 1271 1075 0901 0749 0618 0505 0409 0329 0262 0207 0162 0125 0096 0073 0055 0041 0031 0023 0016 0012 0008 0006 0004 0003
0.05 4801 4404 4013 3632 3264 2912 2578 2266 1977 1711 1469 1251 1056 0885 0735 0606 0495 0401 0322 0256 0202 0158 0122 0094 0071 0054 0040 0030 0022 0016 0011 0008 0006 0004 0003
0.06 4761 4364 3974 3594 3228 2877 2546 2236 1949 1685 1446 1230 1038 0869 0721 0594 0485 0392 0314 0250 0197 0154 0119 0091 0069 0052 0039 0029 0021 0015 0011 0008 0006 0004 0003
0.07 4721 4325 3936 3557 3192 2843 2514 2206 1922 1660 1423 1210 1020 0853 0708 0582 0475 0384 0307 0244 0192 0150 0116 0089 0068 0051 0038 0028 0021 0015 0011 0008 0005 0004 0003
0.08 4681 4286 3897 3520 3156 2810 2483 2177 1894 1635 1401 1190 1003 0838 0694 0571 0465 0375 0301 0239 0188 0146 0113 0087 0066 0049 0037 0027 0020 0014 0010 0007 0005 0004 0003
0.09 4641 4247 3859 3483 3121 2776 2451 2148 1867 1611 1379 1170 0985 0823 0681 0559 0455 0367 0294 0233 0183 0143 0110 0084 0064 0048 0036 0026 0019 0014 0010 0007 0005 0003 0002
81
