Inverso da potência das distâncias

O método do inverso da potência das distâncias é um procedimento para estimação ou interpolação determinista (por oposição a um processo estocástico) em recurso a um conjunto de dados separados, geralmente no espaço. As localizações de valor desconhecido são calculadas em recurso a uma média ponderada pelo inversos das distâncias dessa mesma localização à localização de valores conhecidos (regra geral os dados reais).^[1] É habitualmente utilizado em sistemas de informação geográfica e geoestatística e por esse motivo comparado com outros dos métodos de interpolação utilizados nessa área como é o caso da krigagem ou interpolação por vizinho mais próximo (também chamado, por vezes, método dos polígonos de influência). A sua formulação habitual é utilizando a potência de 2, sendo nesse caso denominado inverso do quadrado da distância.

Nota histórica

No laboratório de Harvard para computação gráfica e análise espacial, começando em 1965, vários cientistas reconsideraram, entre outras coisas, as metodologias utilizadas em sistemas de informação geográfica. Howard Fisher concebeu um novo software de mapeamento chamado SYMAP com o objectivo de melhorar a interpolação. Divulgou o seu trabalho entre os seus colegas dos quais vários vieram a participar nos eventos do laboratório. Um desses colegas, Donald Shepard, decidiu melhorar o procedimento do SYMAP, sendo os seus desenvolvimentos documentados num artigo em 1968.^[2] O algoritmo de Shepard foi influenciado pelas experiências de William Warntz e outros colegas do mesmo laboratório que trabalhavam em análise espacial e conduziam experiências com a exponenciação da distância. Devido a isso Shepard usou o inverso da potência distância como ponderador na interpolação.

Definição

O cálculo do valor u numa dada posição x em recurso às amostras $u_{i}=u(x_{i})$ for $i=0,1,...,N$ utilizando o método do inverso da potência das distâncias é (segundo Shepard, 1968)^[2]:

u(\mathbf {x} )=\sum _{i=0}^{N}{\frac {w_{i}(\mathbf {x} )u_{i}}{\sum _{j=0}^{N}{w_{j}(\mathbf {x} )}}},

onde

w_{i}(\mathbf {x} )={\frac {1}{d(\mathbf {x} ,\mathbf {x} _{i})^{p}}}

Repare-se que no cálculo do $w$ o valor de $p$ é a potência utilizada. Se, a título de exemplo, for substituído por 2 passar-se-ia a designar por método do inverso do quadrado da distância.
Note-se também que o cálculo do ponderador numa localização onde exista uma amostra é indeterminado dado a distância ser zero e consequentemente com ponderador dividido por zero.
Quanto maior a potência mais influência será dada às amostras mais próximas do ponto a estimar (ver figura a seguir).

Interpolação segundo método de Shepard para diferentes potências p, a partir de pontos dispersos numa malha.

Algoritmo do Inverso da potência das distâncias (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 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.

def IPD(dados,blocos,p):
    resultado = np.zeros(blocos)
    for x in xrange(resultado.shape[0]):
        for y in xrange(resultado.shape[1]):
             dist = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             if 0 in dist:
                 ind = np.where(dist == 0)
                 resultado[x,y] = dados[ind[0][0],2]
             else:
                 resultado[x,y] = np.sum((1/dist**p)*dados[:,2])/np.sum(1/dist**p)
    return resultado

Numa versão simplificada poderiamos ter o seguinte procedimento:

def IPD(dados,blocos,p):
    resultado = np.zeros(blocos)
    for x in xrange(resultado.shape[0]):
        for y in xrange(resultado.shape[1]):
             dist = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             resultado[x,y] = np.sum((1/dist**p)*dados[:,2])/np.sum(1/dist**p)
    return resultado

No entanto repare-se que o algoritmo não está preparado para conceber casos onde a distância do bloco da matriz ao ponto é 0 (sendo por isso o inverso da distância de valor indeterminado). Assim inserindo a seguinte condição if/else podemos prever esta situação perguntando se existe algum valor 0 no vector calculado de distâncias entre pontos dos dados e bloco a calcular:

             dist = np.sqrt((x-dados[:,0])**2+(y-dados[:,1])**2)
             if 0 in dist:
                 ind = np.where(dist == 0)
                 resultado[x,y] = dados[ind[0][0],2]
             else:
                 resultado[x,y] = np.sum((1/dist**p)*dados[:,2])/np.sum(1/dist**p)

Se existir tal caso então o algoritmo procura a posição do 0 e devolve esse valor como valor estimado (sendo por isso o valor da própria amostra onde está localizado), caso contrário é feito o procedimento de calculo habitual. Para melhor compreensão do algoritmo descrito aqui é importante que se compreenda o que algumas das funções fazem:

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, no nosso caso em particular onde o vector de distâncias tem o valor de 0. Existe a possibilidade de haver mais do que uma posição onde isto aconteça devolvendo portanto vários índices. Por esse motivo o resultado do np.where(), ind, é chamado como ind[0][0], o primeiro índice do primeiro vector.
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 uma potência de dois resultou, o seguinte por comparação com a imagem original de onde foram retirados os dados:

Imagem objectivo

Interpolação

Interpolação utilizando o inverso do quadrado da distância implementado no código acima.

A utilização de diferente potências (1 a 6) tem uma influência significativa no modelo final que se pode aproximar mais ou menos da imagem objectivo:

Potência 1

Potência 2

Potência 3

Potência 4

Potência 5

Interpolação utilizando o inverso do quadrado da distância implementado no código acima.

Evidentemente no sentido de se atingir uma imagem objectivo (nesta experiência por nós conhecida) não só é importante a representação visual como também a reprodução estatística dos dados. Não existe qualquer ferramenta matemática que force à reprodução do histograma na interpolação por inverso da potência das distâncias. Mesmo na eventualidade disso se concretizar nada garante que os dados usados na interpolação tenham a mesma distribuição univariada (histograma) que a imagem objectivo (ou área de estudo) como podemos ver na comparação seguinte entre a distribuição dos dados e distribuição da imagem objectivo:

Analisando a reprodução do histograma nas interpolações por oposição à imagem objectivo podemos observar que para este caso em particular quanto maior a potência mais aproximado é o histograma:

Potência 1

Potência 2

Potência 3

Potência 4

Potência 5

Interpolação utilizando o inverso do quadrado da distância implementado no código acima.

Isto não implica que a melhor reprodução da imagem objectivo será com potências elevadas mas sim que dependerá do caso de estudo. Entre excepções de notar estão:

A dimensão da área de estudo em relação à amostragem.
A semelhança entre a distribuição estatística (entre outros atributos) dos dados utilizados na interpolação e a realidade (cuja análise dependerá evidentemente da experiência do modelador).

Discussão

Existem várias modificações ao procedimento original de Shepard, geralmente por modificando o método de cálculo dos ponderados ou inserindo a noção de anisotropia na estimação à semelhança da krigagem.

Ver também

Referências

↑ Soares, A. (2006), "Geoestatística para as ciências da Terra e do Ambiente" (2006), Lisboa: Instituto Superior Técnico
↑ ^a ^b Shepard, Donald (1968). «A two-dimensional interpolation function for irregularly-spaced data». Proceedings of the 1968 ACM National Conference. pp. 517–524. doi:10.1145/800186.810616

[1] Soares, A. (2006), "Geoestatística para as ciências da Terra e do Ambiente" (2006), Lisboa: Instituto Superior Técnico

[shepardArticle-2] Shepard, Donald (1968). «A two-dimensional interpolation function for irregularly-spaced data». Proceedings of the 1968 ACM National Conference. pp. 517–524. doi:10.1145/800186.810616

[1]

[2]