Mostrar código-fonte de Krigagem

'''Krigagem''' ou '''krigeagem''' (do [[língua inglesa|inglês]] '''Kriging''') é um método de [[Regressão (estatística)|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|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''' ('''B'''est '''L'''inear '''U'''nbiased '''E'''stimator - 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 ==
[[Ficheiro:variograma-exemplo.png|thumb|300px|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:
:<math>v\,</math>: variância
:<math>c_0\,</math>: ''nugget''
:<math>a\,</math>: efeito pepita
:<math>c_0+c\,</math>: variância assintótica
:<math>h\,</math>: distância de separação
==== Linear ====
Este modelo não apresenta '''silo''' e é muito simples. Sua curva pode ser representada por:
:<math>v= c_0+ch\,</math>

==== Esférico ====
A forma esférica é a mais utilizada e possui silo. Sua forma é definida por:
:<math>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}</math>

==== Exponencial ====
A curva do variograma exponencial respeita a seguinte equação:
:<math>v=c_0+c (1-e^\frac{-h}{b})\,</math>
==== Gaussiano ====
A forma gaussiana é dada por:
:<math>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}</math>

== 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.
:<math>v(h=d_{ip})=\frac{1}{2n}\sum_{i=1}^{n}(f_i-f_p)^2</math>
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:
:<math>F(x,y)=\sum_{i=1}^{n}w_if_i</math>
onde <math>n</math> é o número de amostras obtidas, <math>f_i</math> é o valor obtido no ponto <math>i</math> e <math>w_i</math> é o peso designado ao ponto <math>i</math>.
A fim de obter os pesos de cada um dos <math>n</math> pontos, para cada um deles é realizado um cálculo de <math>w_1, w_2, ..., w_n</math>. Tal procedimento depende do tipo de '''Krigagem''' que está a ser utilizado. Salienta-se a seguinte notação:
:<math>w_j\,</math>: peso do ''j-ésimo'' ponto
:<math>S(d_{ij})\,</math>: valor da semi-variância de <math>d_{ij}</math>
:<math>\lambda\,</math>: 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:
:<math>
\begin{cases}
w_1S(d_{11})+w_2S(d_{12})+...+w_nS(d_{1n})+\lambda=S(d_{1p}) \\
w_1S(d_{21})+w_2S(d_{22})+...+w_nS(d_{2n})+\lambda=S(d_{2p}) \\
\wr \\
w_1S(d_{n1})+w_2S(d_{n2})+...+w_nS(d_{nn})+\lambda=S(d_{np}) \\
w_1+w_2+...+w_n=1
\end{cases}
</math>

==== 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 <math>\lambda</math> 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:
:<math>
\begin{cases}
w_1S(d_{11})+w_2S(d_{12})+...+w_nS(d_{1n})=S(d_{1p}) \\
w_1S(d_{21})+w_2S(d_{22})+...+w_nS(d_{2n})=S(d_{2p}) \\
\wr \\
w_1S(d_{n1})+w_2S(d_{n2})+...+w_nS(d_{nn})=S(d_{np})
\end{cases}
</math>

=== Obtendo o Ponto Interpolado ===
Ao obtermos os valores de <math>w_1, w_2, ..., w_n</math>, calcula-se o valor de <math>f_p</math>:
:<math>f_p=w_1f_1+w_2f_2+...+w_nf_n\,</math>
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]]:
*[http://www.nbb.cornell.edu/neurobio/land/OldStudentProjects/cs490-94to95/clang/kriging.html 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):

<syntaxhighlight lang="Python">
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
</syntaxhighlight>

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).

<syntaxhighlight lang="Python">
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
</syntaxhighlight>

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:

<syntaxhighlight lang="Python">
    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
</syntaxhighlight>

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).

<syntaxhighlight lang="Python">  
    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()
</syntaxhighlight>

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.

<syntaxhighlight lang="Python">  
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
</syntaxhighlight>

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:

<syntaxhighlight lang="Python">  
    K[-1,:] = 1
    K[:,-1] = 1
    K[-1,-1] = 0
</syntaxhighlight>

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:

[[File:Sparse data.png|thumb|400px|center|Dados utilizados no exemplo de estimação.]]

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:

{{multiple image
   | align     = center
   | width     = 400
   | footer    = Interpolação utilizando a krigagem simples e normal implementada no código acima.
   | image1    = Objective image - Channel reservoir.png
   | alt1      = Objectivo
   | caption1  = Imagem objectivo
   | image2    = Simple_kriging_80_50.png
   | alt2      = Previsão
   | caption2  = Interpolação por krigagem simples com variograma (80,50)
   | image3    = Ordinary_kriging_80_50.png
   | alt3      = Previsão
   | caption3  = Interpolação por krigagem normal com variograma (80,50)
  }}

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:

{{multiple image
   | align     = center
   | width     = 400
   | footer    = Interpolação utilizando a krigagem simples com dois modelos de anisotropia.
   | image1    = Simple_kriging_80_50.png
   | alt1      = Previsão
   | caption1  = Interpolação por krigagem simples com variograma (80,50)
   | image2    = Simple_kriging_50_80.png
   | alt2      = Previsão
   | caption2  = Interpolação por krigagem simples com variograma (50,80)
  }}

== Ligações externas ==
*[http://www.dpi.inpe.br/spring/usuario/krigeagem.htm SPRING - Geoestatística - Krigeagem]
*[https://web.archive.org/web/20070504144046/http://gasa.dcea.fct.unl.pt/gasa/tig/AulaT5/Geoestatistica.htm Modelos de Dados Geográficos]
*[https://web.archive.org/web/20100825012551/http://www.ems-i.com/gmshelp/Interpolation/Interpolation_Schemes/Kriging/Kriging.htm Kriging(ems-i)]

{{Portal3|Probabilidade e Estatística}}

[[Categoria:Estatística]]
[[Categoria:Geoestatística]]