Krigagem
# 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):
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).
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:
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).
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.
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:
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:
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: