Krigagem

Krigagem ou krigeagem (do inglês Kriging) é um método de regressão usado em geoestatística para aproximar ou interpolar dados. A teoria de Kriging foi desenvolvida a partir dos trabalhos do seu inventor, Daniel G. Krige, pelo matemático francês Georges Matheron, no começo dos anos sessenta. Na comunidade estatística, também é conhecido como “Processo Gaussiano de Regressão”. A estimação com base em apenas um atributo insere-se no âmbito da Krigagem; a estimação de um atributo à custa de outros atributos insere-se no âmbito da Co-krigagem.

Introdução

Krigagem poderá ser entendido como uma predição linear ou uma forma da Inferência bayesiana. Parte do princípio que pontos próximos no espaço tendem a ter valores mais parecidos do que pontos mais afastados. A técnica de Krigagem assume que os dados recolhidos de uma determinada população se encontram correlacionados no espaço. Isto é, se num aterro de resíduos tóxicos e perigosos a concentração de Zinco num ponto p é x, é muito provável que se encontrem resultados muito próximos de x quanto mais próximos se estiver do ponto p (princípio da geoestatística). Porém, a partir de determinada distância de p, certamente não se encontrarão valores aproximados de x porque a correlação espacial pode deixar de existir.

Considera-se o método de Krigagem do tipo BLUE (Best Linear Unbiased Estimator - Melhor Estimador Linear não-Viciado): é linear porque as suas estimativas são combinações lineares ponderadas dos dados existentes; é não enviesada pois procura que a média dos erros (desvios entre o valor real e o valor estimado) seja nula; é a melhor porque os erros de estimação apresentam uma variância (variância de estimação) mínima. O termo Krigagem abrange um conjunto de métodos, sendo os mais usuais os seguintes:

Tipos de Krigagem

Krigagem Simples

Assume que as médias locais são relativamente constantes e de valor muito semelhante à média da população que é conhecida. A média da população é utilizada para cada estimação local, em conjunto com os pontos vizinhos estabelecidos como necessários para a estimação.

Krigagem Normal

As médias locais não são necessariamente próximas da média da população usando-se apenas os pontos vizinhos para a estimação. É o método mais usado em problemas ambientais.

Co-krigagem

É uma extensão do anterior a situações em que duas ou mais variáveis são espacialmente dependentes e a variável que se quer estimar não está amostrada com a intensidade com que estão as outras variáveis dependentes, utilizando-se os valores destas e as suas dependências para estimar a variável requerida.

Conceitos matemáticos

O Método de Krigagem serve-se de diversas teorias explanadas na estatística. No entanto, para deixarmos mais claras as teorias de estatística usadas e mais direcionadas ao escopo deste texto, explicaremos alguns conceitos.

Semi-variância e semi-variograma

A semi-variância é a medida do grau de dependência espacial entre duas amostras. A magnitude da semi-variância entre dois pontos depende da distância entre eles, implicando em semi-variâncias menores para distâncias menores e semi-variâncias maiores para distâncias maiores. O gráfico das semi-variâncias em função da distância a um ponto é chamado de Semi-variograma. A partir de uma certa distância a semi-variância não mais aumentará com a distância e estabilizar-se-á num valor igual à variância média, dando a esta região o nome de silo ou patamar (sill). A distância entre o início do semi-variograma e o começo do silo recebe o nome de range ou amplitude ou alcance. Ao extrapolarmos a curva do semi-variograma para a distância zero, podemos chegar a um valor não-nulo de semi-variância. Este valor recebe o nome de Efeito Pepita (Nugget Effect).

Modelos de Variograma

No Método de Krigagem normalmente são usados quatro tipos de variogramas. Neles, são usadas as seguintes variáveis:

v\,

: variância

c_{0}\,

: nugget

a\,

: efeito pepita

c_{0}+c\,

: variância assintótica

h\,

: distância de separação

Linear

Este modelo não apresenta silo e é muito simples. Sua curva pode ser representada por:

v=c_{0}+ch\,

Esférico

A forma esférica é a mais utilizada e possui silo. Sua forma é definida por:

v={\begin{cases}c_{0}+c[1.5({\frac {h}{a}})-0.5({\frac {h}{a}})^{3}],&{\mbox{se }}h<a\\c_{0}+c,&{\mbox{se }}h>a\end{cases}}

Exponencial

A curva do variograma exponencial respeita a seguinte equação:

v=c_{0}+c(1-e^{\frac {-h}{b}})\,

Gaussiano

A forma gaussiana é dada por:

v={\begin{cases}c_{0}+c(1-e^{\frac {-h^{2}}{a^{2}}}),&{\mbox{se }}h<a\\c_{0}+c,&{\mbox{se }}h>a\end{cases}}

O Método de Krigagem

Determinação do Semivariograma

Toma-se como base a simulação de um sistema de duas dimensões (2D) que contém um número finito de pontos onde é possível a medição de uma grandeza qualquer. Após a aquisição destes dados, iniciar-se-á a interpolação por Kriging buscando alcançar uma maior resolução. O primeiro passo é construir um semivariograma experimental. Para tal, calcula-se a semivariância de cada ponto em relação aos demais e insere-se no gráfico da semivariância pela distância.

v(h=d_{ip})={\frac {1}{2n}}\sum _{i=1}^{n}(f_{i}-f_{p})^{2}

A partir deste gráfico estima-se o modelo de variograma que melhor se aproxima da curva obtida. O efeito pepita pode estar presente no semivariograma experimental e deve ser considerado. Determinado o modelo do semivariograma a ser usado, inicia-se a fase de cálculos. Sendo o semivariograma uma função que depende da direção, é natural que apresente valores diferentes conforme a direção, recebendo este fenómeno o nome de Anisotropia. Caso o semivariograma apresente uma forma semelhante em todas as direções do espaço, só dependendo de h, diz-se que a estrutura é Isotrópica, i. e., sem direções privilegiadas de variabilidade.

Cálculo dos Pesos

Considere, para o cálculo da Krigagem, a seguinte fórmula:

F(x,y)=\sum _{i=1}^{n}w_{i}f_{i}

onde $n$ é o número de amostras obtidas, $f_{i}$ é o valor obtido no ponto $i$ e $w_{i}$ é o peso designado ao ponto $i$ . A fim de obter os pesos de cada um dos $n$ pontos, para cada um deles é realizado um cálculo de $w_{1},w_{2},...,w_{n}$ . Tal procedimento depende do tipo de Krigagem que está a ser utilizado. Salienta-se a seguinte notação:

w_{j}\,

: peso do j-ésimo ponto

S(d_{ij})\,

: valor da semi-variância de

d_{ij}

\lambda \,

: variável temporária

Krigagem Normal

Neste caso é utilizada a média local dos pontos amostrados. Por conseguinte, deve-se normalizar a média dos pesos. Consequentemente, tem-se um resultado mais preciso do que a Krigagem Simples. Utilizar-se-ão as seguintes equações para a determinação dos valores dos pesos no p-ésimo ponto:

{\begin{cases}w_{1}S(d_{11})+w_{2}S(d_{12})+...+w_{n}S(d_{1n})+\lambda =S(d_{1p})\\w_{1}S(d_{21})+w_{2}S(d_{22})+...+w_{n}S(d_{2n})+\lambda =S(d_{2p})\\\wr \\w_{1}S(d_{n1})+w_{2}S(d_{n2})+...+w_{n}S(d_{nn})+\lambda =S(d_{np})\\w_{1}+w_{2}+...+w_{n}=1\end{cases}}

Krigagem Simples

Para este caso, utiliza-se a média de todos os dados. Implica-se, portanto, em não se normalizar a média local dos pesos, como no caso anterior. Assim, teremos quase que a mesma equação, exceto pela exclusão de $\lambda$ e pela última equação. A característica principal deste método é a geração de gráficos mais lisos e mais esteticamente suaves. Deve-se salientar que este caso é menos preciso que o caso anterior. Os valores dos pesos para o p-ésimo ponto serão dados por:

{\begin{cases}w_{1}S(d_{11})+w_{2}S(d_{12})+...+w_{n}S(d_{1n})=S(d_{1p})\\w_{1}S(d_{21})+w_{2}S(d_{22})+...+w_{n}S(d_{2n})=S(d_{2p})\\\wr \\w_{1}S(d_{n1})+w_{2}S(d_{n2})+...+w_{n}S(d_{nn})=S(d_{np})\end{cases}}

Obtendo o Ponto Interpolado

Ao obtermos os valores de $w_{1},w_{2},...,w_{n}$ , calcula-se o valor de $f_{p}$ :

f_{p}=w_{1}f_{1}+w_{2}f_{2}+...+w_{n}f_{n}\,

Desta maneira, calcula-se o valor interpolado para todos os pontos desejados. Ressalta-se que somente devem ser utilizados os valores adquiridos acima.

Interpolando Outros Pontos

A obtenção do valor interpolado em um outro ponto requer a repetição de todos os cálculos realizados a partir da obtenção do modelo de variograma. Desta forma, para aumentarmos a resolução que é pretendida, deve-se recorrer à métodos matemáticos para a resolução computacional. Diversos códigos foram desenvolvidos para esta resolução, mas um dos melhores algoritmos pode ser obtido no link abaixo. Ele fora desenvolvido inicialmente para a linguagem Fortran, porém foi recodificado para C com a ajuda da biblioteca fortran2c e apresenta-se totalmente em C:

Kriging Interpolation Algorithm in C

Algoritmo da krigagem (Python)

A implementação bidimensional que se segue é feita na linguagem Python (versão 2.7.2) com recurso à biblioteca NumPy tendo por esse motivo o seguinte cabeçalho de importações (está considerado também a biblioteca matplotlib usada para visualização):

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

Repare-se que a primeira linha de importação refere-se à possibilidade de uma divisão entre inteiros poder ser lida como podendo ter um resultado decimal não sendo de maior consequência na transcrição do seguinte procedimento para qualquer outra linguagem (inclusive o próprio Python em versões mais recentes). A função que se segue faz, com efeito, a interpolação por krigagem simples com recurso a uma matriz de dados na qual a primeira coluna é a coluna da posição X, a segunda posição Y, e terceira valor da amostra. O último argumento trata-se de uma tupla ou lista com o valor da amplitude na direção X (posição 0) e direção Y (posição 1).

def KS(dados,blocos,variograma):
    resultado = np.zeros(blocos)
    
    cov_angulos = np.zeros((dados.shape[0],dados.shape[0]))
    cov_distancias = np.zeros((dados.shape[0],dados.shape[0]))
    K = np.zeros((dados.shape[0],dados.shape[0]))
    for i in xrange(dados.shape[0]-1):
        cov_angulos[i,i:]=np.arctan2((dados[i:,1]-dados[i,1]),(dados[i:,0]-dados[i,0]))
        cov_distancias[i,i:]=np.sqrt((dados[i:,0]-dados[i,0])**2+(dados[i:,1]-dados[i,1])**2)
    for i in xrange(dados.shape[0]):
        for j in xrange(dados.shape[0]):
            if cov_distancias[i,j]!=0:
                amp=np.sqrt((variograma[1]*np.cos(cov_angulos[i,j]))**2+(variograma[0]*np.sin(cov_angulos[i,j]))**2)
                K[i,j]=dados[:,2].var()*(1-np.e**(-3*cov_distancias[i,j]/amp))
    K = K + K.T
    
    for x in xrange(resultado.shape[0]):
        for y in xrange(resultado.shape[1]):
             distancias = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             angulos = np.arctan2(y-dados[:,1],x-dados[:,0])
             amplitudes = np.sqrt((variograma[1]*np.cos(angulos[:]))**2+(variograma[0]*np.sin(angulos[:]))**2)
             M = dados[:,2].var()*(1-np.e**(-3*distancias[:]/amplitudes[:]))
             W = LA.solve(K,M)
             resultado[x,y] = np.sum(W*(dados[:,2]-dados[:,2].mean()))+dados[:,2].mean()
    return resultado

O algoritmo acima feito não usa qualquer motor de procura ou limite de amostras. Desta maneira todas as amostras são usadas no sistema de krigagem pelo que apenas é necessário calcular a matriz M, com o valor de variogram entre todas as amostras, uma vez. Isto é feito na seguinte parcela de código:

    cov_angulos = np.zeros((dados.shape[0],dados.shape[0]))
    cov_distancias = np.zeros((dados.shape[0],dados.shape[0]))
    K = np.zeros((dados.shape[0],dados.shape[0]))
    for i in xrange(dados.shape[0]-1):
        cov_angulos[i,i:]=np.arctan2((dados[i:,1]-dados[i,1]),(dados[i:,0]-dados[i,0]))
        cov_distancias[i,i:]=np.sqrt((dados[i:,0]-dados[i,0])**2+(dados[i:,1]-dados[i,1])**2)
    for i in xrange(dados.shape[0]):
        for j in xrange(dados.shape[0]):
            if cov_distancias[i,j]!=0:
                amp=np.sqrt((variograma[1]*np.cos(cov_angulos[i,j]))**2+(variograma[0]*np.sin(cov_angulos[i,j]))**2)
                K[i,j]=dados[:,2].var()*(1-np.e**(-3*cov_distancias[i,j]/amp))
    K = K + K.T

De resto apenas é necessário calcular a matriz M por cada nó da malha de resultados de maneira a resolver o sistema K*W=M (o objectivo é obter W, a matriz de pesos). Após o cálculo dos pesos basta multiplicar os mesmos aos dados e proceder à soma final (com atenção à forma de cálculo em krigagem simples onde: resultado = Soma(W*dados - média) + média). Repare-se que nesta krigagem simples poderia haver um outro argumento que é a média imposta pelo utilizador que posteriormente é utilizado no cálculo do ponto a estimar (nesta implementação é automáticamente considerado a média dos dados de entrada).

  
    for x in xrange(resultado.shape[0]):
        for y in xrange(resultado.shape[1]):
             distancias = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             angulos = np.arctan2(y-dados[:,1],x-dados[:,0])
             amplitudes = np.sqrt((variograma[1]*np.cos(angulos[:]))**2+(variograma[0]*np.sin(angulos[:]))**2)
             M = dados[:,2].var()*(1-np.e**(-3*distancias[:]/amplitudes[:]))
             W = LA.solve(K,M)
             resultado[x,y] = np.sum(W*(dados[:,2]-dados[:,2].mean()))+dados[:,2].mean()

Para fazer krigagem normal as modificações são poucas. A matriz K tem mais uma coluna e mais uma linha e a matriz M tem mais um elemento também. A multiplicação dos pesos (sem considerar o último valor) pelos dados dá-nos o valor estimado em cada bloco.

  
def OK(dados,blocos,variograma):
    resultado = np.zeros(blocos)
    
    cov_angulos = np.zeros((dados.shape[0],dados.shape[0]))
    cov_distancias = np.zeros((dados.shape[0],dados.shape[0]))
    K = np.zeros((dados.shape[0]+1,dados.shape[0]+1))
    for i in xrange(dados.shape[0]-1):
        cov_angulos[i,i:]=np.arctan2((dados[i:,1]-dados[i,1]),(dados[i:,0]-dados[i,0]))
        cov_distancias[i,i:]=np.sqrt((dados[i:,0]-dados[i,0])**2+(dados[i:,1]-dados[i,1])**2)
    for i in xrange(dados.shape[0]):
        for j in xrange(dados.shape[0]):
            if cov_distancias[i,j]!=0:
                amp=np.sqrt((variograma[1]*np.cos(cov_angulos[i,j]))**2+(variograma[0]*np.sin(cov_angulos[i,j]))**2)
                K[i,j]=dados[:,2].var()*(1-np.e**(-3*cov_distancias[i,j]/amp))
    K = K + K.T
    K[-1,:] = 1
    K[:,-1] = 1
    K[-1,-1] = 0
    
    for x in xrange(resultado.shape[0]):
        for y in xrange(resultado.shape[1]):
             distancias = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             angulos = np.arctan2(y-dados[:,1],x-dados[:,0])
             amplitudes = np.sqrt((variograma[1]*np.cos(angulos[:]))**2+(variograma[0]*np.sin(angulos[:]))**2)
             M = np.ones((dados.shape[0]+1,1))
             M[:-1,0] = dados[:,2].var()*(1-np.e**(-3*distancias[:]/amplitudes[:]))
             W = LA.solve(K,M)
             resultado[x,y] = np.sum(W[:-1,0]*(dados[:,2]))
    return resultado

A última linha da matriz K é populada por 1 excepto o último valor (que é zero). O mesmo para a última coluna da mesma matriz como se pode ver nesta parcela de código:

  
    K[-1,:] = 1
    K[:,-1] = 1
    K[-1,-1] = 0

Para melhor compreensão do algoritmo descrito aqui é importante que se compreenda o que algumas das funções fazem:

np.arctan2() - devolve o ângulo em recurso à distância nas coordenadas em Y e X.
np.sum() - devolve a soma de todos os elementos do vector a que se aplica.
np.sqrt() - calcula a raiz quadrada de um dado valor (ver formulação para o cálculo da Distância).

Na utilização da função o utilizador deverá introduzir como argumentos a matriz de dados (onde as colunas devem ser X, Y , valor, por esta ordem: 0,1,2), uma tupla com o número de blocos na dimensão X e Y, por exemplo: (150,200), e finalmente o valor da potência que se pretende utilizar (2 para o inverso do quadrado das distâncias). No exemplo seguinte a estimação vai ser feita numa matriz de 150 blocos na direção X e 200 blocos na direção Y com recurso à função transcrita acima e ao seguinte conjunto de dados:

Da função, utilizando a krigagem simples e normal, resultou o seguinte por comparação com a imagem original de onde foram retirados os dados:

Imagem objectivo

Interpolação por krigagem simples com variograma (80,50)

Interpolação por krigagem normal com variograma (80,50)

Interpolação utilizando a krigagem simples e normal implementada no código acima.

A alteração do modelo de variograma imposto tem influência no resultado final como podemos ver nos dois modelos utilizados na implementação acima da krigagem simples que deram origem aos seguintes resultados:

Interpolação por krigagem simples com variograma (80,50)