Variograma

O variograma ( $2\gamma (x,y)$ ) é uma função que mede a variação do valor de uma variável em relação às restantes da mesma amostragem. Embora seja, de facto, uma derivação da medição de dispersão estatística variância é comumente utilizado em estatística espacial devido a contextualizar esta medição com a dimensão espacial considerando, geralmente mas não obrigatoriamente, a distância entre amostras e/ou a orientação delas.

Definição

Se $\mu =E(Z)$ é o valor esperado (média) de uma dada amostragem $Z$ então a sua dispersão estatística pode ser calculada com a variância:

\operatorname {var} (Z)=\operatorname {E} ((Z-\mu )^{2})

No entanto esta medição não nos dá qualquer informação sobre a sua dispersão espacial pelo que é utilizado o método do variograma para, re-amostrando as populações de acordo com a distância entre pares de pontos, calcular a variância considerando apenas os pares de amostras que se encontram a uma distância $h$ . Assim o variograma é dado pela seguinte fórmula:

2\gamma (x,x+h)={\text{var}}(Z(x)-Z(x+h))=E\left(|(Z(x)-E(x))-(Z(x+h)-E(x+h))|^{2}\right)\quad

Onde o valor de $\gamma (x,x+h)$ é chamado de semi-variograma (o semi-variograma é também muitas vezes, no campo da geoestatística, designado variograma). Se assumirmos que a média é constante por todo o campo espacial considerado então o variograma é:

2\gamma (x,x+h)={\frac {1}{N(h)}}\sum _{N(h)}[Z(x)-Z(x+h)]^{2}\quad

Onde o $N(h)$ é o número de pares de pontos à distância de $h$ . Repare-se que ao fazer este cálculo para re-amostragens de pontos com diferentes distâncias $h$ é possível fazer uma previsão da variância espacial expectável para a amostragem $Z$ . O critério da distância não é o único que pode ser utilizado na caracterização espacial da variância (variograma), também a orientação pode ser usada na re-amostragem da mesma população. Se a orientação de uma direção sobre um dado referencial for dada pelos ângulos " $\alpha$ , $\beta$ " então a re-amostragem para o cálculo do valor do variograma só poderá ter em conta pares de pontos que se encontrem com uma orientação espacial em relação ao referencial usado de " $\alpha$ , $\beta$ " à distància de $h$ . Quando a orientação das re-amostragens para o cálculo do variograma é desconsiderada, este é comumente designado por variograma omni-direcional caso contrário trata-se de um variograma direcional com orientação " $\alpha$ , $\beta$ ".

Variograma experimental

O variograma experimental (ou variograma empírico) é calculado a partir dos valores amostrais da variável consoante a distância e a direcção consideradas, visualizando o resultado por meio de um gráfico de dispersão. Assumindo que re-amostramos uma dada população $Z$ e calculámos o valor do semi-variograma para as distâncias de $h_{1}$ , $h_{2}$ , $h_{3}$ , $h_{4}$ e $h_{5}$ , obtemos a previsão experimental da variação espacial da população $Z$ .

Importa referir que calcular o variograma de todos os pares de amostras a uma distância $h$ ou calcular o valor esperado (média) dos vários variogramas individuais é equivalente:

E[\gamma (h)]={\frac {1}{2|N(h)|}}\sum _{(i,j)\in N(h)}E[|Z(x_{i})-Z(x_{j})|^{2}]={\frac {1}{2|N(h)|}}\sum _{(i,j)\in N(h)}2\gamma (x_{j}-x_{i})

Devido a vários factores (nomeadamente a escassez de dados) muitas vezes se dão tolerância aos critérios de cálculo do variograma (distância e orientação). Assim um par de amostras que se encontre suficientemente perto de ter uma distância $h$ entra para o cálculo do valor do variograma desse mesmo $h$ ( $h\pm \Delta h$ ). Igualmente podemos dizer que um par de amostras que encontre suficiente perto de ter uma orientação " $\alpha$ , $\beta$ " entra para o cálculo do valor do variograma dessa mesma orientação " $\alpha$ , $\beta$ " (" $\alpha \pm \Delta \alpha$ , $\beta \pm \Delta \beta$ "). Na figura abaixo podemos ver como vários pontos próximos das suas distâncias $h$ (a verde) são usados no cálculo do variograma final (a preto).

Outros critérios são usados no cálculo do variograma experimental, dependendo do objectivo ou software utilizado, como por exemplo:

distância mínima considerada (cutoff distance).
distância máxima considerada (cutoff distance).
limite angular (bandwidth).
número de classes (number of lags).
tamanho da classe (lag distance).
tolerância ao tamanho da classe (lag tolerance).

Em termos gerais podemos dizer que o variograma experimental corresponde ao variograma real calculado a partir dos dados implicando, necessariamente, tratar-se de uma análise discreta dos dados pelo que não é possível saber o valor de variograma para todas as distâncias $h$ .

Modelos de variograma

Em alguns campos, como é o caso da geoestatística (geralmente ao utilizar-se métodos como, por exemplo, a krigagem) é necessário conseguir-se calcular o valor do variograma para qualquer distância o que leva a ter de se ajustar um modelo matemático aos dados experimentais. Isto é feito em recurso ao um patamar imposto pelo utilizador (que, regra geral, é igual à variância da amostragem). Este patamar aparece como uma linha no variograma experimental e irá fazer com que modelo ajustado não exceda o limite imposto convergindo para o mesmo:

Embora muitos modelos matemáticos possam ser ajustados a um variograma experimental os mais usados são (Soares, 2006)^[1]:

Modelo exponencial:

\gamma (h)=C_{0}+C_{1}[1-e^{(-3h/a)}]\qquad

Modelo esférico:

\gamma (h)=\left\{{\begin{matrix}C_{0}+C_{1}[1.5{\frac {h}{a}}-0.5({\frac {h}{a}})^{3}],&para\quad h\leq a\\C,&para\quad h>a\end{matrix}}\right.

Modelo gaussiano:

\gamma (h)=C_{0}+C_{1}[1-e^{(-3h^{2}/a^{2})}]\qquad

0 ajuste do modelo é feito com base em três parâmetros:

Efeito pepita (nugget effect, $C_{0}$ ).
Patamar (sill, $C=C_{0}+C_{1}$ ).
Amplitude (range, $a$ ).

Em teoria o valor do variograma é nulo, $\gamma (h)=0$ , quando $h=0$ , no entanto na prática existe uma diferença no valor para o mais pequeno $h$ pelo qual possa ser quantificado $\gamma (h)$ . Quando este valor é elevado então assume-se existir uma grande variabilidade à pequena escala implicando que $\gamma (h)$ não tende para zero quando $h$ tende para zero. São nestes casos que importa ajustar o modelo ao variograma experimental considerando o efeito a pequenas escalas dado pela constante $C_{0}$ designada por efeito pepita. Na figura seguinte pode-se ver um modelo esférico ( ${\color {Green}verde}$ ), exponencial ( ${\color {Blue}azul}$ ) e gaussiano ( ${\color {Purple}purpura}$ ) ajustados a um variograma experimental com a mesma amplitude e efeito pepita nulo (origem dos modelos estão no ponto "0,0").

Algoritmo do variograma (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

De notar que o código seguinte introduzido em três funções não segue todas os parâmetros geralmente utilizados pelos softwares de geoestatística. A título de exemplo num variograma a direção de 90º é igual à direção de -90º o que não está previsto nesta implementação. Serve apenas como demonstração do cálculo de uma variograma experimental e posterior adequação de um modelo.

def variograma_experimental(dados):
    angulos = np.zeros((dados.shape[0],dados.shape[0]))
    distancias = np.zeros((dados.shape[0],dados.shape[0]))
    angulos[:,:]=np.NAN
    distancias[:,:]=np.NAN
    semivariancias = np.zeros((dados.shape[0],dados.shape[0]))
    for i in xrange(dados.shape[0]-1):
        angulos[i,i:]=np.arctan2((dados[i:,1]-dados[i,1]),(dados[i:,0]-dados[i,0]))
        distancias[i,i:]=np.sqrt((dados[i:,0]-dados[i,0])**2+(dados[i:,1]-dados[i,1])**2)
        semivariancias[i,i:]=((dados[i:,2]-dados[i,2])**2)/2
    angulos = (angulos*180)/np.pi
    for i in xrange(angulos.shape[0]):
        for j in xrange(angulos.shape[1]):
            if angulos[i,j]<0:
                angulos[i,j] = angulos[i,j] +180
    return angulos, distancias,semivariancias
    
def variograma_direcional(angulos,distancias,semivariancias,direcao,tolerancia,classes):
    ind = np.where((angulos > direcao - tolerancia) & (angulos < direcao + tolerancia))
    distancias_direcionais = distancias[ind]
    semivariancias_direcionais = semivariancias[ind]
    hist = np.histogram(distancias_direcionais,classes)
    resultado = np.zeros((classes-1,2))
    for i in xrange(classes-1):
        ind2 = np.where((distancias_direcionais > hist[1][i]) & (distancias_direcionais < hist[1][i+1]))
        resultado[i,0]=hist[1][i]+(hist[1][1]-hist[1][0])/2
        resultado[i,1]=np.mean(semivariancias_direcionais[ind2])
    return resultado
    
def modelo_variograma(resultado,amplitude,patamar):
    plt.scatter(resultado[:,0],resultado[:,1],color='green',s=130)
    vector_patamar=np.zeros(resultado.shape[0])
    vector_patamar[:]=patamar
    vector_distancia=resultado[:,0].copy()
    vector_distancia[0]=0
    vector_distancia[-1]=vector_distancia[-1]+(vector_distancia[-1]-vector_distancia[-2])
    plt.plot(vector_distancia,vector_patamar,color='red',linewidth=3)
    modelo = patamar*(1-np.e**(-3*np.linspace(0,vector_distancia[-1],1000)/amplitude))
    plt.plot(np.linspace(0,vector_distancia[-1],1000),modelo,color='blue',linewidth=2)
    plt.ylabel('Variograma')
    plt.xlabel('Distancia')
    plt.grid()
    plt.xlim(0,vector_distancia[-1])
    plt.ylim(0,vector_patamar[-1]+0.2*vector_patamar[-1])
    plt.show()

A primeira função cria matrizes com o valor de variograma, ângulos e distâncias entre todos os pontos. Isto é feito em recurso às funções cuja explicação se segue:

np.zeros() - Cria uma matriz com a dimensão introduzida de colunas e linhas preenchida com zeros.
np.NAN - Trata-se de uma valor sem dados, implicando que não deve ser utilizado, apenas preenche espaço necessário.
np.arctan2() - Devolve o ângulo em radianos (é feita posteriormente a conversão para graus) em recurso a coordenada cartesianas.

Após o cálculo do valor do variograma para todos os pares de pontos possíveis é feita selecção em recurso à segunda função cujos argumentos de entrada são as matrizes resultado da função anterior, a direção, tolerância e o número de classes que o variograma deverá ter. Funções importantes a ter em conta para compreensão do algoritmo são:

np.where() - devolve um vector (ou vector de vectores, dependendo da dimensão do objecto) com as posições onde a condição que foi imposta é realizada.
np.histogram() - devolve uma tupla com dois vectores sendo o segundo a separação dos dados de input em cada classe (limite superior e inferior) e o primeiro o número de elementos que está em cada classe. Neste caso em particular apenas o segundo vector é utilizado.
np.mean() - Cálculo da média do vector de input.

A terceira função recebe o resultado da segunda com os valores de variograma e distância escolhidas e elabora um gráfico onde aparece, além destes pontos, a linha do patamar e um modelo exponencial ajustado com efeito pepita zero. Tem como argumentos para além do variograma experimental a amplitude do modelo e o valor do patamar. De notar as funções:

plt.scatter() - Cria um gráfico de pontos.
plt.plot() - Cria um gráfico de linhas.
plt.xlabel() e plt.ylabel() - Insere o nome de cada eixo (X e Y).
plt.xlim() - Dá os limites da janela de visualização para o eixo do X.
plt.grid() - Insere um quadriculado de fundo para o gráfico.
plt.show() - Faz a visualização do gráfico numa janela (até esta instrução o gráfico está invisível).

Usando o seguinte conjunto de dados como exemplo:

Foram calculados os variogramas na direções 0º e 90º que devolveu o seguinte resultado:

Direcção 0º

Direcção 90º

Cálculo do variograma experimental e respectiva adequação do modelo.

Discussão

Em geoestatística são usadas habitualmente três funções para estudar a variabilidade espacial da amostragem que são: covariância, correlograma, e semi-variograma (comumente designado variograma). A figura seguinte mostra o variograma experimental, covariância e correlograma para o mesmo conjunto de dados:

Existe também variogramas cruzados para utilização com métodos de co-estimação como por exemplo a co-krigagem e variogramas de indicatriz para métodos de estimação com variáveis de indicatriz, como por exemplo krigagem da indicatriz.