FILTRAGEM DE NUVEM LASER PARA GERAÇÃO DE MDT POR KRIGAGEM

Laser point cloud filtering for DTM generation using kriging

Resumos

Esse artigo apresenta um método de filtragem de nuvem de pontos LASER para obtenção de um Modelo Digital de Terreno (MDT). O processo de filtragem é realizado com base em uma superfície aproximada gerada a partir de pontos amostrados sobre vias urbanas. Esses pontos são determinados por meio de linhas detectadas na imagem de intensidade do pulso laser via detector de Steger. A principal suposição do método é que o terreno tem um comportamento suave no interior das quadras e, dessa forma, os pontos amostrados ao longo das vias permitem, utilizando o método de interpolação por krigagem, uma representação adequada do terreno no interior das quadras, ou seja, relativamente próxima dos pontos LASER de terreno dessas regiões. Assim, a filtragem é realizada pela verificação da proximidade dos pontos da nuvem LASER original com a superfície aproximada gerada. Por fim, um MDT é obtido da nova amostra pelo método de interpolação por krigagem, melhorando a descrição da superfície. A partir dos experimentos realizados foi possível verificar a viabilidade do método proposto, com resultados de boa coerência visual e indicadores numéricos satisfatórios

Palavras-chave:
MDT; Filtragem; Krigagem

This paper presents a method of filtering point clouds generated by laser scanning, to obtain a Digital Terrain Model (DTM). The filtering process is performed based on an approximated surface obtained from urban road points. These points are sampled in straight lines detected by Steger in the intensity image of the laser pulse. The main assumption of the method is that the ground has smooth behavior inside the block, so the sample laser points collected along the urban roads allow, using the kriging interpolation method, a suitable representation of the land inside the block, that is, relatively close to the ground point laser in these regions. Thus, filtering is performed by proximity of the original cloud laser points with approximate surface. For thus a DTM is obtained from the new sample by kriging interpolation method, increasing the description of the surface. From the experiments it was possible to verify the feasibility of the proposed method, with results of good visual consistency and satisfactory numerical indicators.

Keywords:
DTM; Filtering; Kriging

1. Introdução

A geração de MDT (Modelo Digital de Terreno), a partir de dados de Varredura a LASER (Light Detection and Ranging) Aerotransportado - VLA vem se tornando comum em razão do desenvolvimento das técnicas de aquisição de dados 3D por VLA e processamento destes dados. O MDT pode ser obtido da nuvem de pontos LASER por meio de um processo de filtragem, que consiste na remoção dos pontos que não pertencem à superfície do terreno, ou seja, aqueles que se encontram acima do solo, tais como construções, árvores, entre outros.

Uma das dificuldades na geração de MDT a partir de dados de VLA é a variação na densidade de pontos obtida do solo, uma vez que pode haver regiões com abundância de pontos enquanto outras com escassez em virtude de vegetação e construções que impedem a incidência do LASER no solo. Com isto torna-se necessário um processo cuidadoso de filtragem da nuvem de pontos, de forma a garantir que os pontos retidos sejam de terreno.

Trabalhos como de Assunção (2010), destacam as diferentes técnicas de filtragem e classificação para geração automática de MDT. Um ponto em comum nestas técnicas é a dificuldade em identificar pontos que possam gerar uma superfície inicial para o processo de filtragem.

O método proposto no presente trabalho faz uso de ferramentas de PDI (Processamento Digital de Imagens) e Geoestatística para a filtragem da nuvem de pontos LASER, a fim de produzir um MDT acurado. As ferramentas de PDI são usadas para selecionar pontos sobre vias de tráfego de veículos na imagem de intensidade do pulso laser. Esses pontos são amostras para gerar uma superfície aproximada por krigagem, a qual é usada no processo de filtragem. O principal pressuposto do método para a aquisição dessa superfície aproximada é o comportamento suave do terreno no interior das quadras, que permite utilizar os pontos de vias urbanas como amostras e o interpolador por krigagem para obtenção de uma representação adequada do terreno no interior das quadras. A suavidade do terreno no interior das quadras é esperada para o sucesso do método e está fundamentada na expectativa de baixa aclividade/declividade dos terrenos urbanos, conforme Lei Federal N° 6766 de 19/12/79, que dispõe sobre o Parcelamento do Solo Urbano. Essa superfície aproximada gerada é chamada nesse trabalho como MDT inicial adotado como parâmetro de filtragem da nuvem de pontos original, para a obtenção do MDT final. A partir dos experimentos realizados foi possível verificar a viabilidade do método proposto.

1.1 Modelo Digital de Elevação

O Modelo Digital de Elevação (MDE) é a forma mais utilizada para representar uma superfície de maneira digital com base em um conjunto de pontos com coordenadas tridimensionais. Assim, um MDE consiste em uma representação discreta de uma superfície topográfica, na qual as elevações do terreno podem ser representadas computacionalmente por um conjunto de pontos regularmente distribuídos (Wolf e Dewwitt, 2000).

Um MDE pode se subdividir em Modelo Digital do Terreno (MDT) ou Modelo Digital de Superfície (MDS). O MDT consiste em um MDE no qual as informações de pontos de elevação se restringem a superfície topográfica, compreendendo apenas o terreno e corpos d’agua. E o MDS refere-se a um MDE que apresenta informações de elevação de pontos na superfície terrestre e de objetos que estão sobre ela (tais como árvores, edificações, etc) (El-Sheimy et al., 2005). A obtenção de um MDE pode ser realizada por meio de técnicas fotogramétricas, como também, através de levantamento por sistemas GNSS (Global Navigation Satellite System), estação total, entre outros. Essas técnicas consistem da aquisição de uma malha de pontos com coordenadas de terreno que permitam a modelagem almejada (El-Sheimy et. al, 2005). Uma alternativa que se tornou viável atualmente é a aquisição de dados por meio de sistemas de VLA. Para a elaboração de modelos de superfícies existem vários processos. O Triangulated Irregular Network (TIN) e a malha regular são os métodos mais usados para representar superfícies.

Quando se trata da obtenção de uma malha regular a escolha da função de interpolação é de fundamental importância para se obter uma boa precisão do modelo. Uma função interpoladora deve possibilitar a reprodução de uma superfície contínua e um tempo computacional não restritivo, além de possuir propriedades matemáticas de interesse para uma dada aplicação.

a) Método de Interpolação por Krigagem

A interpolação por krigagem é um método geoestatístico. A geoestatística tem por base a teoria das variáveis regionalizadas. O comportamento dessas variáveis é baseado em alguns pressupostos, como o de estacionariedade, que propõe que o fenômeno é descrito como homogêneo na região em que se pretende fazer estimativas; o de ergodicidade que enuncia que a esperança referente a média de todas as possíveis realizações da variável é igual a média de uma única realização dentro de um certo domínio; e finalmente o pressuposto da hipótese intrínseca que diz que as diferenças entre valores apresentam fraco incremento, ou seja, são localmente estacionárias (Landim et al., 2002).

O método de interpolação por krigagem é fundamentado na teoria das variáveis regionalizadas e, dessa forma, leva em consideração tanto a distância entre as amostras como o seu agrupamento. Pela técnica krigagem estima-se valores médios e uma medida de acuracidade dessa estimativa. Os pesos a serem associados às amostras são calculados com base na distância entre a amostra e o ponto estimado, na continuidade espacial e no arranjo geométrico do conjunto (Bettini, 2007). A determinação desses pesos está diretamente relacionada à modelagem da continuidade espacial obtida a partir do semivariograma.

O semivariograma pode ser obtido experimentalmente, com uso do esquema de amostragem em duas dimensões apresentado na Figura 1, em que Z(x ₁) e Z(x ₁+ h) correspondem ao valor em uma posição cujos componentes são (x ₁,y ₁) e (x₂ ,y₂ ) , respectivamente; h é um vetor que aponta de (x ₁,y ₁) para (x₂ ,y₂ ), com modulo igual à distância entre esses pontos (Camargo et al., 2002).

Figura 1:
Semivariograma. (a) Gráfico de distribuição espacial, e (b) componentes do semivariograma. Fonte: Camargo et al., 2002

Assim, para determinar experimentalmente o semivariograma é considerado para cada valor de h, todos os pares de pontos Z (x_i ) e Z (x_i +h), separados pela distância ||h|| e a função semivariograma (γ(h)) é calculada a partir da Equação 1 (Camargo et al., 2002).

(1)

com n(h) definido pelo número de pares de observações.

A função γ(h) descreve o desvio esperado nos valores dos pontos Z (x_i ) em função da distância entre pares de pontos. A Figura 1b apresenta o semivariograma experimental, em função de γ(h) e h, e relaciona-o a um comportamento idealizado composto pelos seguintes parâmetros (Camargo et al., 2002):

Alcance (a): representa a distância máxima em que as amostras se apresentam correlacionadas espacialmente;

Patamar: valor de γ(h) correspondente ao seu Alcance (a); e

Efeito pepita (C₀): representa a descontinuidade do semivariograma para distâncias menores que a menor distância entre as amostras.

Dentre várias técnicas de krigagem, a ordinária consiste em uma forma de estimação linear para uma variável regionalizada que atende à hipótese intrínseca, ou seja, não requer o conhecimento prévio da média e assume-se a hipótese de estacionaridade local (Yamamoto e Landim, 2013). A Equação 2 fornece o estimador da krigagem ordinária.

(2)

em que Z(x_i ) são os dados experimentais; λ_i, (i = 1,2,..,n) são os pesos atribuídos a cada valor amostral e n, o número total de dados.

1.2 Filtragem de dados VLA

Um sistema de VLA produz um conjunto de pontos 3D irregularmente espaçados. Esse conjunto de pontos é denominado nuvem de pontos. As técnicas de filtragem de dados de VLA têm por finalidade detectar e eliminar os pontos da nuvem LASER que não pertencem ao terreno. Como consequência, pontos associados com vegetação e outros objetos elevados devem ser eliminados nesse processo. Os processos de filtragem já desenvolvidos, segundo Sithole e Vosselman (2004), podem ser conceituados de acordo com vários critérios, dentre os quais:

1) A estrutura dos dados utilizada: os algoritmos de filtragem utilizam uma nuvem de pontos original (Axelsson, 1999; Sithole, 2001) ou uma malha regular de pontos (reamostrada) (Brovelli, 2002; Elmqvist, 2001);

2) Critérios de definição de vizinhança: os algoritmos operam em uma vizinhança local e na classificação um ou mais pontos são classificados por vez. Essa classificação envolve três técnicas distintas: classificação ponto contra ponto, classificação ponto contra pontos e classificação pontos contra pontos (Sithole e Vosselman, 2004);

3) Medida de descontinuidade: para a separação do terreno de objetos elevados, os algoritmos de filtragem se baseiam em uma medida de descontinuidade. Comumente, são utilizadas, por exemplo, diferença de altura, declividade, menor distância a faces de um TIN, e menor distância a superfícies parametrizadas (Sithole e Vosselman, 2004);

4) Princípios de filtragem:

a) Baseado na declividade: Nesses algoritmos (Sithole, 2001), a declividade entre dois pontos é medida. Se o desnível supera certo limiar, o ponto mais alto é assumido como pertencente a um objeto elevado;

b) Baseado no bloco mínimo: A função discriminante é um plano horizontal delimitando uma região em 3D, na qual os pontos do terreno devem ser encontrados (Wack e Wimmer, 2002);

c) Baseado na superfície parametrizada: semelhante ao anterior, com a diferença de que o plano é substituído por uma superfície parametrizada (Elmqvist, 2001) e Axelsson, 1999); e

d) Baseado na segmentação: esse método (por exemplo, Brovelli, 2002) considera que pontos que se agrupam e estão acima de uma vizinhança, pertencem a um objeto. Para tanto, os agrupamentos devem delinear objetos e não faces de objetos.

5) Mecanismo de controle: o iterativo ou o não iterativo (Sithole e Vosselman, 2004); e Natureza da filtragem: trata da remoção ou da recolocação dos pontos filtrados no conjunto de dados (Sithole e Vosselman, 2004).

2. Método

O método utilizado nesse trabalho pode ser dividido em três etapas. A primeira está relacionada à obtenção de um MDT inicial por meio da interpolação por krigagem em pontos amostrados nas vias urbanas. A segunda etapa consiste na filtragem da nuvem LASER original, com base no MDT obtido na etapa anterior, para gerar a partir dos pontos filtrados um MDT final por krigagem. Na terceira etapa são propostas formas de validação dos resultados através de uma avaliação qualitativa (análise visual), e quantitativa (índices numéricos). As três etapas serão descritas a seguir.

2.1 Geração do MDT Inicial

A nuvem de pontos LASER, que é o dado de entrada, possui pontos com coordenadas tridimensionais (X, Y e Z), no sistema geodésico WGS (World Geodetic System) 84, e também a resposta radiométrica (I) do pulso LASER (Jensen, 2009).

O primeiro passo para a tarefa de obtenção do MDT inicial consiste em adquirir a amostra de pontos pertencentes ao terreno. Para tanto, optou-se por adquirir pontos de vias. Essa tarefa tem por base a imagem de intensidade, produzida a partir das coordenadas X, Y e I dos pontos LASER, e as características das vias nessas imagens. O asfalto possui índice de reflexão baixo na faixa do infravermelho (faixa de operação da maioria dos sistemas VLA), e menor que o índice de reflexão dos outros objetos presentes nas cenas urbanas (Wehr e Lohr, 1999). Assim as vias aparecem escuras nas imagens de intensidade, correspondendo no histograma da imagem ao pico mais a esquerda. Dessa forma, para o processo de isolamento das regiões de vias basta aplicar uma limiarização, com base em um limiar adquirido através da análise do histograma da imagem de intensidade. O resultado da limiarização consiste em uma imagem em que as regiões de vias correspondem aos pixels pretos e os demais objetos passam a ser codificados na cor branca. Em seguida, na imagem limiarizada é aplicado um filtro de mediana, para a eliminação de ruídos provocados por pequenos agrupamentos isolados de pontos presentes na nuvem LASER, destoando da feição de interesse.

O próximo passo consiste em aplicar o detector de linhas de Steger (Steger,1996) na imagem resultante para extrair segmentos de vias. Esse processo é seguido da eliminação de retas com comprimento mínimo, para minimizar os falsos positivos (segmentos presentes no resultado que não representam vias).

Assim, como a imagem de intensidade está registrada com a nuvem de pontos LASER, basta separar os pontos da nuvem que pertencem aos segmentos de vias detectados nessa imagem.

Em seguida, as retas obtidas são percorridas ponto a ponto em sua extensão, comparando a diferença em altura com base em um limiar d_alt . Os pontos fora desse limiar são considerados não pertencentes às vias e desprezados, eliminando descontinuidades em altura com relação aos pontos de vias.

Por fim, o MDT inicial é obtido, em malha regular, com uso do método de interpolação por krigagem. Nesse trabalho utiliza-se a interpolação por krigagem uma vez que se assume que o comportamento do terreno no interior das quadras é suave e esse método de interpolação conserva melhor essa característica, por manter a continuidade do terreno (Botelho, 2007).

2.2 Filtragem da nuvem LASER pelo método de superfície e geração do MDT

O processo de filtragem tem por base o MDT inicial obtido na etapa anterior e o conceito de célula da malha regular, que é a região delimitada por quatro pontos vizinhos. A busca por pontos de terreno na nuvem LASER é realizada célula a célula deste MDT através da aproximação do segmento de superfície associado a cada célula por duas faces triangulares adjacentes. Assim, o método consiste em inicialmente, encontrar os pontos que pertencem à (célula), com base nas coordenadas planimétricas da célula e do ponto. Em seguida, determinar a qual das faces triangulares da célula pertence cada ponto. Esse procedimento é realizado com base no cálculo do posicionamento relativo (Figura 2), entre o ponto Ps e o segmento de reta (definido pela diagonal formada pelo ponto superior esquerdo Pse e o ponto inferior direito Pid da célula). Assim, conforme pode ser observado na Figura 2, para determinar se o ponto está a esquerda ou a direita do segmento de reta, basta calcular a área formada pelo triângulo PsePsPid pela Equação 3. Se esta for positiva o ponto Ps está a esquerda, negativa o ponto Ps está a direita e nula os pontos Pse, Ps e Pid estão alinhados. (Davis Jr e Queiroz, 2005).

(3)

Definida a face triangular da célula a qual pertence o ponto selecionado, determina-se, para cada ponto (Ps (Xs, Ys, Zs)), o ponto Pt(Xs, Ys, Zt = f(Xs, Ys)) contido na face triangular t que tem como suporte a equação do plano f , a distância dst entre Ps e Pt. Para o caso dos pontos estarem alinhados (Figura 2c) o valor de Zt é obtido através de uma interpolação linear com base nos valores de Z dos pontos extremos da reta.

Figura 2:
Posicionamento relativo de ponto e segmento de reta, com a) ponto a esquerda; b) ponto a direita e c) ponto alinhado com o segmento de reta

Assim, são eliminados os pontos da nuvem LASER que contêm distâncias dst superiores a um limiar pré-definido. Por fim, um MDT final é obtido por krigagem utilizando os pontos filtrados e tidos como de terreno.

2.3 Análise dos Resultados

Para a avaliação dos resultados foram adotados critérios qualitativos e quantitativos. O critério qualitativo se pautou na interpretação visual dos resultados por parte de um operador humano, que observou a suavidade do modelo digital de terreno gerado e o seu poder de descrição do relevo. Já para a análise quantitativa foram calculados os seguintes indicadores numéricos: Média (Equação 4), Desvio-Padrão (Equação 5) e Raiz Quadrada do Erro Quadrático Médio - RMSE, sigla em inglês (Equação 6), por meio das discrepâncias altimétricas (ΔZ). As discrepâncias altimétricas foram obtidas pela diferença entre os pontos de referência, colhidos manualmente na nuvem de pontos LASER (com o auxílio da imagem de intensidade), e os correspondentes no MDT final. Os pontos de referência foram selecionados no meio das quadras, ou seja, afastados das vias onde foram retiradas as amostras iniciais.

(4)

(5)

(6)

Com $\stackrel{—}{\text{Δ}Z}$ sendo a média das discrepâncias altimétricas (ΔZ_i ) e n é o tamanho da amostra (número de pontos de controle).

Ao final, foi realizada uma avaliação de tendência espacial da distribuição das discrepâncias altimétricas com base numa análise de superfícies de tendência (Landim, 1997). Essa avaliação é realizada pela aplicação de análise de variância nos resíduos da representação dessa distribuição espacial como sendo um plano médio com outro que representa uma superfície com tendência espacial. Permite verificar que o procedimento não foi espacialmente tendencioso e, portanto, as discrepâncias distribuem-se aleatoriamente no espaço.

3. Resultados Experimentais

No desenvolvimento das diferentes etapas do método foram utilizados os softwares Visual Studio C++, MVTec Halcon 12.0 e Surfer 8.0. Para a experimentação foram utilizadas duas regiões da área urbana da cidade de Curitiba/PR, denominadas Área 1 e Área 2, com áreas de 444.158 m² e 489.240 m², respectivamente. A nuvem LASER utilizada foi adquirida por um sistema aerotransportado e possui densidade aproximada de 2 pontos/m². A Figura 3a e a Figura 3b apresenta as correspondentes imagens de intensidade geradas com uma resolução espacial de 0,5m.

Figura 3:
Imagem de intensidade a) Área 1 e b) Área 2

Na Figura 4 podem ser observadas as imagens de intensidade limiarizadas. O limiar utilizado foi o valor 25 de tom de cinza em uma escala de 0 (preto) a 255 (branco), obtido com base na análise do histograma destas imagens. Verifica-se o destaque das vias em razão de sua resposta espectral.

Figura 4:
Imagem de intensidade limiarizada a) Área 1 e b) Área 2

Apesar das vias estarem bem evidenciadas ainda ocorrem ruídos que afetarão o processo de amostragem. Dessa forma, foi aplicada a suavização por filtro de mediana na imagem de intensidade limiarizada, com elemento estruturante retangular de tamanho 3. O resultado dessa filtragem pode ser observado na Figura 5.

Figura 5:
Imagem de intensidade limiarizada suavizada a) Área 1 e b) Área 2

Na imagem de intensidade limiarizada e suavizada foi aplicado o detector de Steger, seguido da filtragem de retas com o comprimento de até 10 pixels. Após a eliminação de retas pelo comprimento mínimo, as coordenadas de terreno (X, Y, Z) dos pontos que formam cada uma das retas remanescentes foram recuperadas, possibilitando a eliminação de descontinuidades em altura, conforme previsto no método, com base no limiar d_alt = 0,5m, resultando nas Figuras 6 (a) e (b).

Figura 6:
Localização das amostras a) Área 1 e b) Área 2

Por descreverem vias os pontos amostrados estão alinhados e agrupados e, dessa forma, não seguem um padrão ideal de distribuição, o que é evidenciado pela ausência de amostras no interior das quadras. No entanto, assumindo que há um comportamento suave no interior das quadras, tais pontos são representativos do terreno.

Na Figura 7, podem ser observados os variogramas experimentais para as amostras. O modelo que melhor se ajustou ao variograma, para os dois casos, foi o de “potência”, o que indica que os dados apresentam uma tendência, pois não possui um patamar (Yamamoto e Landim, 2013). Assim, aplicou-se o método de krigagem universal para modelar a tendência e a krigagem ordinária nos resíduos (Landim, 1997 - p. 183; Andriotti, 2003 - p. 149; Deutsch and Journel, 1998 - p. 69). Para realizar essa tarefa é necessário preliminarmente modelar uma superfície de tendência e utilizar os resíduos no processo de inferência por krigagem. Para tanto, ajustou-se uma superfície de tendência de 1º grau por MMQ e o respectivo mapa de resíduos (Figura 8).

Figura 7:
Variograma das amostras de vias da (a) Área 1 e (b) Área 2

Figura 8:
Mapa da superfície de tendência de 1° grau das Áreas 1 e 2 (a, b) e resíduos (c, d).

A partir da Figura 8, pode-se observar que os resíduos apresentam valores baixos, o que demonstra que o terreno apresenta um comportamento de relevo semelhante próximo ao da superfície ajustada pela função polinomial de 1° grau.

Obtidos os resíduos, foram gerados os variogramas de superfície para as amostras 1 e 2 (Figuras 9 (a) e (b)), para verificar se os dados possuem comportamento anisotrópico ou isotrópico.

Figura 9:
Variograma de superfícies das amostras a) 1 e b) 2.

A partir da observação da distribuição do variograma de superfície, para o primeiro caso (Figura 9a) verificou-se o comportamento isotrópico e para o segundo caso (Figura 9b) um comportamento anisotrópico, que foi simplificado para isotrópico, pois os modelos são similares e, dessa forma, restringiu-se o raio de busca na interpolação. Os variogramas apresentados na Figura 10 foram gerados com 25 intervalos com incrementos de 14,8 e 16,8, respectivamente para as amostras 1 e 2. Verificou-se que o modelo que melhor se ajustou foi o de “potência”, indicando tendência nos resíduos. Esse comportamento pode ser resultado de características inerentes do tipo de dado utilizado (relevo), que em se tratando de pequenas extensões urbanas apresenta variações suaves e contínuas. E, assim, sendo a tendência uma característica decorrente da superfície analisada foi adotado o modelo “potência” no ajuste dos resíduos, com parâmetros de inclinação de 0,0897 e efeito pepita de 0,761 para a amostra 1 e inclinação de 0,0344 para a amostra 2.

Figura 10:
Variograma dos resíduos e respectivo modelo da (a) Amostra 1 e (b) Amostra 2

Definido o modelo de função aplicou-se o processo de krigagem ordinária nos resíduos (Figura 11).

Figura 11:
Superfície por krigagem dos resíduos da (a) Área 1 e (b) Área 2

Ao final, a superfície obtida pela krigagem ordinária dos resíduos foi adicionada à superfície de tendência como correção, produzindo os MDTs iniciais para as áreas 1 (Figura 12a) e 2 (Figura 12b).

Figura 12:
MDT Inicial da (a) Área 1 e (b) Área 2

Os MDTs iniciais obtidos foram então utilizados no processo de filtragem da nuvem de pontos LASER, promovendo assim o adensamento das amostras (Figura 13). Nessa etapa foi utilizado o limiar dst = 0,5m, previsto no método.

Figura 13:
Amostra de pontos adensada da (a) Área 1 e (b) Área 2

Conforme esperado pelo método, obteve-se um adensamento dos pontos descritores do terreno, com novos pontos agrupados nas vias e nas quadras.

De posse dessa nova amostra, o processo de interpolação foi repetido. Por se tratar das mesmas áreas o comportamento de tendência foi novamente verificado e tratado da mesma forma, obtendo-se os MDTs finais para as áreas 1 (Figura 14a) e 2 (Figura 14b).

Figura 14:
MDT Final da (a) Área 1 e (b) Área 2

Seguindo a forma de avaliação descrita na subseção 3.3, considerando o critério qualitativo, o MDT resultante tanto para a Área 1, quanto para a Área 2, foi capaz de descrever o relevo de forma bastante satisfatória, mantendo a suavidade e continuidade esperada para o relevo urbano, sendo possível observar o traçado das vias. Para a análise numérica, a partir de pontos de referência previamente selecionados foi gerado os valores constantes na Tabela 1 para as Áreas 1 e 2. Foram utilizados 23 pontos para a Área 1 e 32 pontos para a Área 2, distribuídos de maneira homogênea nos dois locais.

Tabela 1:
Resumo de índices estatísticos dos MDT.

Para a Área 1 o RMSE do MDT final foi de 0,216 m, 34,35% menor que o RMSE do MDT Inicial (0,329 m), enquanto que para a Área 2 o RMSE do MDT Final ficou em 0,169 m, 44,23% menor que o do MDT Inicial (0,303m). A redução no valor do RMSE para ambas as áreas era esperada em razão do método proposto, que promove um adensamento de pontos da amostra. A Área 2 apresentou uma redução maior (44,23%), em razão de seu relevo ser bastante suave.

Por fim foram gerados gráficos (Figura 15), das discrepâncias entre os pontos de controle e seus homólogos no MDT Inicial e Final, que ajudam na interpretação dos resultados obtidos pelo método.

Figura 15:
Discrepâncias do MDT Inicial e Final da (a) Área 1 e (b) Área 2.

Na Figura 15, a representação visual permite observar que houve melhoras significativas nos pontos com as maiores discrepâncias, sendo que na Área 1 um ponto de controle que apresentava um erro de 0,865 m no MDT Inicial foi reduzido para 0,570 m no MDT Final, enquanto na Área 2 outro ponto com erro de 0,795 m no MDT Inicial foi reduzido para 0,464 m no Final.

Por fim, foi realizada a verificação do ajuste da superfície de tendência às discrepâncias encontradas. Assim, para ambas as áreas foram determinados os parâmetros de uma superfície polinomial de 1°grau (aplicação do método dos mínimos quadrados - MMQ) e em seguida foi feita a validação do modelo através da análise de variância. Nessa verificação, foi realizado um teste de hipótese com nível de significância de (() 5% em que é feita uma análise para verificar se ocorre ajuste significativo da superfície aos dados (H₀) ou não (H₁). Para tanto, foi criada uma tabela de análise de variância (ANOVA - ANalysis Of VAriance). Mais detalhes podem ser encontrados em Landim (1997).

Tabela 2:
Análise de variância para verificação do ajuste de superfície.

O valor obtido para F_c (F calculado) foi comparado ao F_t (F tabelado), encontrado na tabela de distribuição F, de modo que, se F_c fosse maior que F_t a hipótese nula é aceita, ao contrário, é rejeitada. O valor tabelado para a Área 1 é igual a F_t = 3,493 e Área 2 igual a F_t = 3,328 , assim, a hipótese nula foi rejeitada e constatou-se que para ambos os casos não há tendência na distribuição dos erros.

4. Conclusão

Nesse trabalho foi desenvolvido um método para filtragem de nuvem de pontos LASER para a obtenção de um MDT. A filtragem é baseada em uma superfície aproximada obtida com uso de pontos de vias extraídos na imagem de intensidade e o método de interpolação krigagem. Assim, o processo de filtragem é realizado a partir da proximidade dos pontos da nuvem LASER com a superfície aproximada. Um MDT final é obtido pela nova amostra, resultante da filtragem, pelo método de interpolação krigagem, melhorando a descrição da superfície.

Para avaliar o método proposto foram realizados experimentos com nuvens de pontos de duas áreas. Os resultados alcançados indicaram a eficiência do método proposto, em relação aos aspectos tratados a seguir. Os pontos amostrados inicialmente foram obtidos sobre vias, garantindo uma boa confiabilidade na componente altimétrica, porém provocando uma amostragem irregular e bastante agrupada. Para contornar este problema a escolha da interpolação por Krigagem se mostrou acertada, pois a krigagem leva em consideração tanto a distância entre as amostras quanto o seu agrupamento (clustering). A análise (qualitativa e quantitativa), dos resultados indicaram uma melhora significativa entre o MDT Inicial e o MDT Final, resultado da utilização do MDT inicial para um adensamento controlado das amostras, via filtragem, característica do método proposto.

O diferencial desta proposta em relação ao que existe na literatura refere-se à forma como é gerada a aproximação inicial da superfície do terreno, utilizando pontos sobre vias, cujas caracteristicas de extensão, largura, distribuição em ambiente urbano e reposta radiométrica registradas pela varredura LASER (X, Y, Z, I), garantem uma região potencial para coleta de amostras.

Para trabalhos futuros pretende-se aplicar o método em áreas maiores, a fim de verificar se com o aumento da variabilidade do terreno (em cenário urbano), o comportamento de tendência se mantém. Também pretende-se variar a densidade da nuvem de pontos LASER

AGRADECIMENTOS

Os autores agradecem às agências FAPESP(Fundação de Amparo à Pesquisa do Estado de São Paulo) (Processo 2012-22332-2) e FAPEMAT (Fundação de Amparo à Pesquisa do Estado de Mato Grosso) pelo apoio financeiro executado por meio de bolsas de doutorado

REFERÊNCIAS BIBLIOGRÁFICAS

• Andriotti, J. L. S. Fundamentos de Estatística e Geoestatística. São Leopoldo: Editora Unisinos, 2003. 165 p.
• Assunção, M. G. T. Desenvolvimento de Métodos de Filtragem e Classificação de Pontos LASER para a Geração Automática do Modelo Digital do Terreno. Dissertação de mestrado. Universidade Federal de Pernambuco, 2010.
• Axelsson, P. “Processing of LASER Scanner Data - Algorithms and Applications.” ISPRS Journal of Photogrammetry and Remote Sensing, 54 (1999): 138- 47.
• Bettini, C. “Conceitos básicos de Geoestatística.” Em Geomática: modelos e aplicações ambientais. Editado por Meirelles, M. S. P.; Câmara, G.; Almeida, C.M. 193 - 234. Brasília: Embrapa, 2007.
• Botelho, M. F. Modelagem tridimensional de edificações usando dados do sistema laser scanner e imagem orbital de alta. Tese de Doutorado, Universidade Federal do Paraná, 2007.
• Brovelli, M. A., Cannata, M., Longoni, U. M. Managing and processing LASER data within GRASS. Paper presented at GRASS Users Conference, Trento, Italy, September 11-13, 2002.
• Camargo, E. C. G.; Fucks, S. D.; Câmara, G. “A Análise Espacial de Superfícies.” Em Análise Especial de Dados Geográficos. Editado por Druck, S.; Carvalho, M.S.; Câmara, G. e Monteiro, A. M. V. 79-108. Brasília: Embrapa, 2002.
• Davis Jr, C. A., Queiroz, G. R. Algoritmos Geométricos e Relacionamentos Topológicos. In: Casanova, M. A., Câmara, G., Davis Jr., C. A., Vinhas, L., Queiroz, G. R. (Eds.)Bancos de Dados Geográficos Curitiba (PR): EspaçoGeo, 2005, p. 53-92.
• Deutsch, C. V.; Journel, A. G. GSLIB: Geoestatistical Software Library and User's Guide. New York, Oxford University Press, 1998. 369p.
• El-Sheimy, N.; Valeo, C.; Habib, A. Digital terrain modeling: Acquisition, manipulation and applications. Calgary: Artech House, 2005.
• Elmqvist, N. Ground Estimation of Lasar Radar Data Using Active Shape Models. Paper presented at OEEPE Workshop on Airborne LASER scanning and Interferometric SAR for Detailed Digital elevation Models, Frankfurt am Main, Germany, March 1-3, 2001.
• Jensen, J. R. Sensoriamento Remoto do Ambiente: Uma Perspectiva em Recursos Terrestres. Tradução da 2a. edição. Epiphanio, J. C. N. (org.). São José dos Campos: Parêntese Editora, 2009, 672 p.
• Landim, P. M. B. Análise estatística de dados geológicos. Editora: UNESP, 1997.
• Landim, P. M. B.; Sturaro, J. R.; Monteiro, Rubens C. Krigagem ordinária para situações com tendência regionalizada. Rio Claro: DGA, IGCE, UNESP, 2002.
• Sithole, G. Filtering of LASER Altimetry Data Using a Slope Adaptive Filter. Paper presented at XXXIV International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Annapolis, October 22-24, 2001.
• Sithole, G., Vosselman, G. “Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds.” ISPRS Journal of Photogrammetry and Remote Sensing , 59 (2004): 85-101.
• Steger, C. Extracting Lines Using Differential Geometry and Gaussian Smoothing. Paper presented at International Archives of Photogrammetry and Remote Sensing, Vienna, Austria, 1996 , vol. 31, Part B3, pp. 821-826.
• Wack, R., Wimmer, A. Digital Terrain Models from Airborne LASER Scanner Data-a Grid Based Approach. Paper presented at International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Graz, Austria, September 9-13, 2002.
• Wehr, A.; Lohr, U. Airborne laserscanning-an introduction and overview. ISPRS Journal of Photogrammetry and Remote Sensing . v. 54, n. 2-3, pp.68-82, 1999.
• Wolf, P. R.; Dewitt, B. A. Elements of Photogrammetry: with applications in GIS. Boston: McGraw-Hill, 2000.
• Yamamoto, J. K.; Landim, P. M. B. Geoestatística: conceitos e aplicações. Editora: Oficina de Textos, 2013

Datas de Publicação

• Publicação nesta coleção
  Mar 2017

Histórico

• Recebido
  28 Jul 2015
• Aceito
  25 Maio 2016