SciELO - Scientific Electronic Library Online

 
vol.21 issue4Analysis of the Gravimetric Earth Tide at the San Juan Station (Argentina)The photogrammetric spatial resection using quaternions author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Boletim de Ciências Geodésicas

Print version ISSN 1413-4853On-line version ISSN 1982-2170

Bol. Ciênc. Geod. vol.21 no.4 Curitiba Oct./Dec. 2015

http://dx.doi.org/10.1590/S1982-21702015000400043 

Artigos

AVALIAÇÃO AUTOMÁTICA DA ACURÁCIA RELATIVA DE DADOS LIDAR AEROTRANSPORTADO

Automatic assessment of relative accuracy of airborne LiDAR

Gabriel Dresch1  2 

Daniel Rodrigues dos Santos2 

1Diretoria de Serviço Geográfico Exército Brasileiro 1- Brasília - DF - Brasil.Email: dresch@dsg.eb.mil.br

2Programa de Pós-Graduação em Ciências Geodésicas Universidade Federal do Paraná 2 - Curitiba - PR - Brasil. Email: danielsantos@ufpr.br

RESUMO

Neste trabalho é proposto um método automático para o controle da acurácia relativa de dados LiDAR (Light Detection and Ranging) aerotransportados usando superfícies planas, linhas retas e pontos. Primeiramente, é realizada uma filtragem dos dados LiDAR, através do filtro morfológico progressivo e o filtro RANSAC (Random Sample Consensus). Em seguida, as superfícies planas são extraídas automaticamente sobre os telhados encontrados nas faixas adjacentes através do algoritmos crescimento de região. Posteriormente, as correspondências entre as feições extraídas em cada faixa são estabelecidas e os parâmetros de transformação são estimados através de uma variante do método ICP. Finalmente, para avaliar simultaneamente a acurácia plani-altimétrica dos dados LiDAR são empregadas linhas retas extraídas a partir da interseção de planos adjacentes. Para avaliar o método proposto foi conduzido um experimento com dados LiDAR obtidos numa região urbanizada recoberta por 2 faixas com densidade de 1 ponto/m². A avaliação da acurácia relativa desses dados é realizada automaticamente e os resultados obtidos com o método proposto mostraram que a média das distâncias ponto-a-plano é de -0,0019 cm e a raiz quadrada do erro médio quadrático é de 6,04 cm.

Palavras-chave: LiDAR; acurácia relativa; superfícies planas; linhas retas; ICP

ABSTRACT

In this paper an automatic method for assessment of relative accuracy of airborne LiDAR (Light Detection and Ranging) data using surface planes, straight lines and points is presented. Firstly, the LiDAR data is filtered using the progressive morphological filter and RANSAC (Random Sample Consensus) filter. Then, the roof planes are automatically extracted on both strips by region growing algorithm. After, the feature correspondences are established and the transformation parameters are estimated using a variant of ICP algorithm. Finally, in order to evaluate simultaneously planimetric and altimetric relative accuracy, straight lines obtained by intersection of two adjacent surface planes are used. An experiment was conducted to evaluate the proposed method using two overlapping LiDAR strips, whose density is 1 point / m². The assessment of relative accuracy of these data was performed automatically and the results obtained by proposed method, showing values equal -0,0019 cm to mean of the point-to-plane distance and a rmse around 6,04 cm

Keywords: LiDAR; relative accuracy; Surface planes; Straight lines; ICP

1. Introdução

De acordo com Wehr e Lohr (1999) o LiDAR (Light Detection and Ranging) é um sistema de sensoriamento remoto ativo capaz de determinar coordenadas tridimensionais (3D) de pontos sobre a superfície terrestre de forma rápida, precisa e segura por meio da geração, emissão e captura de pulsos LASER (Light Amplification by Stimulated Emission of Radiation) integrados com dados do sistema de posicionamento GNSS (Global Navigation Satellite System) e do sistema de navegação IMU (Inertial Measurement Unit). Por isso, atualmente, esse sistema é empregado em diversos projetos de Engenharia, Geologia e Arquitetura e Urbanismo, tais como, planejamento, hidrografia, mapeamento florestal, monitoramento de áreas costeiras, modelagem 3D de ambientes urbanos, estradas e rodagem, entre outros.

Geralmente, o recobrimento de uma área da superfície física pelo sistema LiDAR é feito por diversas faixas de voo. Os dados são independentemente georreferenciados pela integração dos sensores GNSS/IMU e LASER. Neste procedimento, há um acúmulo de erros sistemáticos e aleatórios nos dados 3D introduzidos pelo mecanismo de varredura do sistema LiDAR e pelas incertezas dos parâmetros de integração das componentes dos sensores supracitados obtidos no processo de calibração, resultando em desalinhamentos angulares e deslocamentos entre as nuvens de pontos 3D relativas a cada faixa de voo. A quantificação e remoção desses erros é um procedimento utilizado pelos produtores de dados LiDAR conhecido como ajustamento de faixas. Além disso, a identificação do deslocamento linear e do desalinhamento angular constitui uma importante parte do processo de avaliação da acurácia feita pelo usuário final. Deste modo, a avaliação da qualidade dos dados LiDAR aerotransportado é de interesse tanto do produtor quanto do usuário. Por isso, é essencial o desenvolvimento de métodos rápidos e confiáveis para o ajustamento de faixas, sendo que tais métodos podem também ser empregados na avaliação da acurácia relativa dos dados LiDAR.

A análise da qualidade de dados LiDAR pode ser relativa ou absoluta. O controle de qualidade que emprega os próprios dados a serem analisados no processo é denominado de avaliação da acurácia relativa, enquanto que na avaliação da acurácia absoluta são utilizados dados externos independentes (López, 2013), tais como, os pontos de controle. Sendo assim, dois aspectos devem ser tratados na avaliação da acurácia relativa de dados LiDAR, isto é: os tipos de erros envolvidos; e o grau de automação a ser empregado no processo.

No primeiro aspecto, como descrito anteriormente, a qualidade dos dados provenientes do LiDAR depende principalmente dos erros aleatórios e sistemáticos oriundos das observações de posição e de orientação obtidas na integração GNSS/IMU, bem como das medidas angulares dos espelhos e das distâncias provenientes do dispositivo de varredura LASER. Já no segundo aspecto, o alto grau de automação também é um fator importante para a avaliação da acurácia relativa, uma vez que implica na redução de tempo de execução e dos custos do projeto. Vale ressaltar que o grande volume de dados gerados pelo LiDAR torna o processo de automação uma tarefa desafiadora e que merece atenção da comunidade científica.

O presente trabalho propõe um método automático para o controle da acurácia relativa de dados LiDAR aerotransportado combinando feições pontuais, lineares e superfícies planas. O método proposto permite avaliar a qualidade relativa da nuvem de pontos LiDAR através da comparação de feições correspondentes extraídas em faixas adjacentes. Neste trabalho também é proposto um método de filtragem de dados LiDAR que permite a aplicação de uma nova variante do algoritmo ICP (Iterative Closest Point) na avaliação da qualidade.

Este trabalho está organizado em quatro seções principais. A seção 2 apresenta os trabalhos relacionados com o tema de estudo. Na seção 3 o método proposto é cuidadosamente apresentado. Os experimentos conduzidos e as discussões dos resultados são apresentados na seção 4. Finalmente, a seção 5 apresenta as conclusões e recomendações para trabalhos futuros.

2. Trabalhos Relacionados

Na última década, diversas abordagens foram propostas para avaliar a acurácia relativa dos dados LiDAR, tendo como foco principal evitar o emprego de alvos pré-sinalizados e de levantamentos em campo. Em geral, as abordagens propostas consistiam em avaliar a acurácia relativa dos dados LiDAR comparando feições correspondentes extraídas automaticamente em faixas LiDAR adjacentes. Killian et al. (1996) adaptaram a técnica de fototriangulação de imagens para concatenar faixas LiDAR adjacentes. A desvantagem da abordagem proposta pelos autores está no fato de que a identificação de pontos homólogos numa nuvem de pontos 3D não apresenta a mesma confiabilidade obtida na identificação de pontos em imagens digitais. Crombaghs et al. (2000) propuseram um método de ajustamento de faixas LiDAR para reduzir discrepâncias altimétricas entre faixas adjacentes. Apesar de eficiente, o procedimento não contempla o tratamento dos erros planimétricos, consideravelmente superiores aos erros altimétricos. Pfeifer et al. (2005) propuseram um método de ajustamento de faixas LiDAR baseado na detecção e estabelecimento de superfícies planas horizontais correspondentes. Os autores empregaram superfícies planas correspondentes para estimar os erros altimétricos, embora nenhuma atenção tenha sido dispensada aos erros planimétricos. Friess (2006) empregou planos correspondentes detectados em faixas adjacentes para estimar correções nas observações oriundas do sistema LiDAR. Segundo os autores, os parâmetros do plano são estimados previamente e mantidos fixos durante o processo de ajustamento de observações. Skaloud e Lichti (2006) estimaram simultaneamente os parâmetros dos planos e os parâmetros derivados da calibração do sistema através de um ajustamento combinado não linear. Já Rentsch e Krzystek (2009) apresentaram um método de ajustamento de faixas com 5 parâmetros que utiliza pontos 2D e 3D extraídos a partir da interseção de linhas de cumeadas. Os parâmetros de transformação são obtidos através de um ajustamento combinado empregando pontos de controle e de enlace, extraídos a partir da interseção de linhas de cumeada.

No estado da arte também são encontrados métodos para avaliação da acurácia absoluta de dados LiDAR usando diferentes tipos de primitivas. Csanyi e Toth (2007) empregaram alvos específicos em voos pré-sinalizados estimando com sucesso a acurácia planimétrica e altimétrica de dados LiDAR. Neste trabalho, os autores ressaltaram que mesmo uma cuidadosa calibração do sistema LiDAR não permite eliminar totalmente os erros sistemáticos. De forma geral, grande parte desses erros sistemáticos pode ser corrigida através do ajustamento de faixas, restando apenas os erros oriundos do sistema de posicionamento GNSS, que não podem ser totalmente eliminados sem a avaliação da acurácia absoluta. Já Toth et al. (2008) afirmaram que embora a solução de Csanyi e Toth (2007) apresente resultados de alta qualidade, o método proposto é limitado por fatores econômicos, já que a confecção, locação e rastreio das coordenadas dos alvos pré-sinalizados demanda tempo e recursos financeiros. Para evitar este problema, os autores sugeriram como controle de campo o emprego de marcas do pavimento de rodovias. Apesar de poderem ser facilmente medidas através de levantamentos geodésicos, o quesito custo ainda esta presente no processo de avaliação, sendo necessária a determinação de dados externos com altíssima precisão.

Alguns autores buscaram alternativas baseadas na avaliação da acurácia relativa, através do emprego de feições correspondentes extraídas em faixas adjacentes. O conceito fundamental desta abordagem reside no fato de que, na ausência de erros, elementos conjugados devem coincidir perfeitamente. Vosselman (2008) analisou a acurácia planimétrica de dados LiDAR através de linhas de cumeada extraídas por interseção de planos obtidos pela transformação de Hough 3D combinado com o algoritmo de crescimento de região. Cada linha extraída fornece a informação da componente do deslocamento perpendicular à linha reta. Neste caso, quando a interseção de planos gera uma linha horizontal, é possível separar as componentes altimétrica e planimétrica do deslocamento entre faixas. Para estimar o erro planimétrico, a posição do centro de cada linha horizontal de uma faixa é comparada em relação à sua linha homóloga na faixa adjacente. Essas linhas provêm da interseção de planos devidamente ajustados, sendo rigidamente reconstruídas. Este método foi automatizado e aplicado por Vosselman (2012) numa região de 100 mil hectares com cerca de 20 bilhões de pontos.

Habib et al. (2010) testaram quatro métodos para avaliação da acurácia relativa de faixas adjacentes: o método ICP, o método ICPatch (Iterative Closest Patch), um método baseado em feições lineares conjugadas e um método baseado em feições planares conjugadas. Os autores concluíram que com exceção do ICP, todos os métodos proporcionam uma estimativa adequada dos erros sistemáticos e aleatórios. Por fim, os autores recomendaram o uso do ICPatch, pois esse procedimento é rápido, não é sensível à densidade de pontos, sendo o que menos depende das características do espaço objeto e exige mínima interação com o usuário. Cabe ressaltar que estas características fazem do algoritmo ICPatch uma alternativa para avaliação de dados extraídos a partir de áreas não urbanizadas. Por outro lado, Sande et al.(2010) afirmaram que a desvantagem do ICPatch é que nem sempre uma TIN é adequada para representar superfícies reais. Já no trabalho de Sande et al.(2010), foram empregados planos extraídos na região de sobreposição entre faixas para estimar erros sistemáticos e aleatórios. Tais planos foram obtidos pelo ajustamento robusto de planos empregando o algoritmo RANSAC (Random Sample Consensus). O modelo apresentado incorpora um grande número de observações confiáveis de distância entre ponto e plano numa estimação por mínimos quadrados. Caso as faixas não possam ser ajustadas apenas por uma translação 3D, uma transformação afim é empregada na estimativa dos parâmetros de transformação. A vantagem deste método é que os parâmetros podem ser obtidos através de um ajustamento paramétrico linear, ou seja, são determinadas de uma única vez, sem necessidade de iterações e linearização do modelo. Outra vantagem do método supracitado é que não há necessidade de valores iniciais aproximados para obtenção da solução.

3. Método

O método proposto neste trabalho está dividido em 3 etapas. A primeira etapa consiste no emprego de dois filtros, isto é: o filtro morfológico progressivo e o filtro RANSAC. A segunda etapa consiste na extração de feições e o estabelecimento de suas correspondências nas nuvens de pontos 3D. Por fim, a avaliação da acurácia relativa dos dados LiDAR é realizada através de uma variante do algoritmo ICP, sendo calculados os parâmetros de transformação entre as faixas e determinadas as precisões altimétrica e planimétrica usando linhas retas.

3.1 Filtro morfológico progressivo e filtragem RANSAC

Primeiramente, o filtro morfológico progressivo proposto por Zhang et al. (2003) é empregado para descartar os pontos que estão no nível do solo, permanecendo apenas pontos pertencentes às vegetações ou edificações. Para eliminar a vegetação remanescente resultante da aplicação do filtro morfológico progressivo é empregado o filtro RANSAC. O RANSAC foi desenvolvido por Fischler e Bolles (1981) e consiste num método para estimação de parâmetros de um modelo a partir de um conjunto de dados que apresenta outliers. O filtro RANSAC foi empregado para remover pontos que não pertencem a regiões planas.

O número t de conjuntos mínimos necessários para extrair um ponto pertencente a uma região plana com 99% de confiança é calculado a partir da seguinte equação:

sendo t o número de conjuntos mínimos, K o tamanho da vizinhança e M a quantidade mínima de pontos necessária para que a região seja considerada plana.

Recomenda-se descartar planos horizontais e verticais para diminuir o número de falsas correspondências e aumentar a eficiência da extração automática de linhas retas.

3.2 Extração de feições e estabelecimento de suas correspondências

Esta etapa consiste em extrair automaticamente superfícies planas e linhas retas nos dados LiDAR. A Figura 1 ilustra o esquema da extração automática dessas feições.

Figura 1: Esquema da extração de feições. 

Na Figura 1, FX1** e FX2** correspondem aos dados de cada faixa pré-processados e filtrados. PL1 e PL2 correspondem aos planos extraídos em cada faixa LiDAR. PT1 e PT2 correspondem aos pontos de cada uma das faixas associados a PL1 e PL2. Cada ponto de PT1 / PT2 está associado a um plano de PL1 / PL2 respectivamente, enquanto LI1 e LI2 correspondem às linhas retas horizontais extraídas a partir da interseção dos planos.

Primeiramente, é realizada a extração de planos empregando o algoritmo de crescimento de regiões proposto por Rabbani et al. (2006). Após a aplicação desse algoritmo, é realizado um processo de refinamento dos dados para eliminar possíveis outliers. A Figura 2 mostra o processo de refinamento dos planos usando o método proposto para filtragem dos dados.

Figura 2: Fluxograma do refinamento de planos. 

De acordo com a Figura 2, os parâmetros dos planos são calculados a priori. A normal do plano é determinada por análise de componentes principais e, posteriormente, normalizada. Já a distância do plano à origem é dada pela média do produto interno entre a normal do plano com cada um dos pontos que o compõem. Em seguida, é calculada a distância de cada ponto ao plano. Pontos que apresentam uma distância superior a ε (limiar pré-estabelecido) são eliminados. Após a eliminação dos pontos, os parâmetros dos planos são recalculados. Esse processo é realizado iterativamente até que as distâncias de todos os pontos ao plano estejam dentro da tolerância estabelecida. Em seguida, é verificado se o plano apresenta um número de pontos superior ao tamanho mínimo do conjunto (M). Se o número de pontos do plano não for suficiente, o mesmo é descartado do processamento.

Em seguida é aplicado o algoritmo Outlier Removal proposto por Rusu (2009) para remover pontos isolados. Para realizar essa tarefa, são estabelecidos dois parâmetros: o raio da esfera d e o número mínimo de vizinhos que devem estar contidos dentro da esfera. Para cada ponto é feita uma busca de vizinhos dentro do raio estabelecido. Se o número de vizinhos for inferior ao valor estabelecido, o ponto é classificado como outlier e é removido do processamento. Em seguida, é novamente verificado se o plano apresenta um número suficiente de pontos, sendo o mesmo descartado caso essa condição não seja atendida. Após a extração dos planos de ambas as faixas, o processo de correspondência é realizado através dos seguintes passos:

  1. É calculada a distância entre o centroide do plano de referência e o centroide de cada plano da faixa de pesquisa. Os planos que apresentarem distância inferior ao valor estabelecido são considerados possíveis candidatos à correspondência;

  2. O plano candidato é descartado se a sua normal (nx, ny, nz) apresentar um ângulo em relação a normal do plano de referência superior a uma tolerância estabelecida;

  3. Se apenas um plano atender essas condições, o mesmo é definido como plano correspondente. Caso haja mais de um plano, o que apresentar menor distância entre os centroides será escolhido. Se nenhum plano atender essas condições, não existe correspondência entre o plano de referência e os de pesquisa;

  4. Por fim, é verificado se há correspondência múltipla, ou seja, se algum plano da faixa LiDAR de pesquisa está relacionado com mais de um plano da faixa de referência. Para evitar ambiguidades e a necessidade de intervenção do usuário, nesses casos, todas as correspondências são eliminadas.

Em seguida, é realizada a extração das linhas retas, sendo dividida em três etapas. A primeira etapa é realizada iterativamente. Em cada iteração, um plano é tomado como plano de referência. Para cada um desses planos, o algoritmo busca por planos vizinhos cuja interseção com o plano de referência determine uma linha horizontal. Para isto, devem ser calculadas as distâncias do centroide do plano de referência em relação aos demais. Os planos considerados vizinhos são aqueles que apresentam uma distância do plano de referência inferior a um valor estabelecido. Sendo assim, é calculado o produto vetorial entre as normais dos planos vizinhos e verificado se a componente vz do produto vetorial apresenta um ângulo inferior a θ em relação à horizontal, como segue:

sendo o vetor diretor da reta, a normal do plano de referência, a normal do plano vizinho e θ a tolerância angular.

Se o vetor resultante do produto vetorial entre as normais apresentar uma inclinação inferior à tolerância pré-estabelecida, a interseção dos planos é considerada horizontal e o número de cada plano é armazenado.

Na segunda etapa, é realizado o cálculo de linhas retas através da interseção de planos adjacentes. Uma linha pode ser representada pela equação paramétrica de uma reta em , a saber:

sendo p um ponto qualquer pertencente à reta , t o parâmetro da reta e o vetor diretor da reta.

A direção da reta é definida pelo produto vetorial entre as normais de dois planos adjacentes. Já o ponto p pode ser qualquer ponto pertencente à reta. Visando obter uma solução unívoca, é aplicada à restrição de que tal ponto deve ser o mais próximo a um determinado ponto de referência. As coordenadas deste ponto são calculadas através dos multiplicadores de Lagrange. A Figura 3 ilustra a representação de uma linha de cumeada definida pela interseção de dois planos (azul e vermelho).

Figura 3: Representação vetorial de uma linha de cumeada. (a) vista em perspectiva; (b) vista superior. 

Na Figura 3, p0 corresponde ao ponto calculado pelo método de Lagrange e corresponde ao vetor diretor da reta . Para fazer a comparação da altitude entre linhas retas correspondentes foi definido um ponto pA (centro aproximado da reta) pertencente à reta que abrange o segmento definido pela linha de cumeada. A posição de pA é definida no espaço pela interseção no plano do segmento definido pelos pontos p1 (centroide do plano vermelho) e p2 (centroide do plano azul) com a projeção da reta no plano xy. O parâmetro t correspondente à posição do ponto pA calculado no plano () e é aplicado na reta em . A coordenada Z do ponto pA é adotada como altitude da linha de cumeada na avaliação da precisão altimétrica.

Por fim, é realizado um refinamento através da eliminação de algumas linhas retas. Esse refinamento tem a finalidade de eliminar linhas imprecisas que possam ter sido selecionadas na primeira etapa. Para isso, cada linha reta é comparada com a sua correspondente e eliminada se não atender os seguintes critérios: a distância entre os centros aproximados deve estar dentro de uma tolerância estabelecida e o ângulo entre essas linhas deve ser inferior à tolerância angular estabelecida.

3.3 Avaliação da Acurácia Relativa

De acordo com Monico et al. (2009) o conceito de acurácia envolve os erros sistemáticos e aleatórios, enquanto o conceito de precisão envolve apenas erros aleatórios. A avaliação da acurácia pode ser feita através do cálculo da raiz quadrado do erro médio quadrático (REMQ), enquanto que para avaliar a precisão é necessário analisar o desvio-padrão das medidas. A Figura 4 ilustra o esquema da avaliação da acurácia relativa.

Figura 4: Esquema da avaliação da acurácia relativa. 

Na Figura 4, PT1 corresponde aos pontos da faixa LiDAR de referência, sendo que cada ponto está associado a um plano da faixa LiDAR de pesquisa (PL2). Os parâmetros de transformação (3 rotações e 3 translações) são calculados pelo algoritmo ICP modificado, proposto neste trabalho, que minimiza o quadrado da distância entre pontos e seus respectivos planos. Em seguida é possível avaliar as distâncias ponto-a-plano antes e após a aplicação dos parâmetros calculados. Já LI1 e LI2 correspondem às linhas retas extraídas de cada uma das faixas. Por serem horizontais, essas linhas permitem a separação das componentes planimétrica e altimétrica, possibilitando a avaliação independente de cada componente. Para isso, os parâmetros calculados são aplicados nas linhas da faixa LiDAR de pesquisa (LI2) e, por fim, é realizada a avaliação da acurácia relativa altimétrica e planimétrica da nuvem de pontos.

3.3.1 Determinação dos Parâmetros de Transformação

Esta etapa consiste em utilizar o algoritmo ICP para encontrar os parâmetros de transformação de modo análogo ao que foi realizado por Sande et al. (2010). Na sua versão original, os autores apresentaram um método que relaciona um conjunto de planos de uma faixa de referência com os pontos correspondentes na faixa de pesquisa. Os parâmetros da transformação entre as duas faixas adjacentes são estimados através do modelo paramétrico do método dos mínimos quadrados (MMQ), o qual minimiza as distâncias entre os pontos e seus planos correspondentes. O método proposto, no presente trabalho, assemelha-se ao método proposto por Sande et al. (2010) nos seguintes aspectos:

(a) Emprega a mesma função erro, ou seja, minimiza o quadrado da distância entre pontos e planos correspondentes;

(b) Utiliza como dados de entrada os planos extraídos de uma faixa e os pontos correspondentes a esses planos da faixa adjacente.

Contudo, apresenta as seguintes diferenças:

(a) Calcula a rotação entre as faixas em vez de estimar uma transformação linear, proporcionando a estimativa direta dos ângulos de desalinhamento do sistema, sem a influência dos demais parâmetros da transformação afim (fatores de escala e fatores de não ortogonalidade) evitando, também, a superestimação do conjunto de parâmetros. A rotação é calculada pelo método proposto por Horn (1987), o qual emprega quatérnios;

(b) Utiliza uma variação do método ICP em vez de utilizar o modelo paramétrico do MMQ.

O algoritmo original do ICP (Besl e Mckay, 1992) minimiza o quadrado da distância entre pontos pseudo-conjugados. Já no método proposto neste trabalho, busca-se minimizar o quadrado da distância entre um conjunto de pontos e planos correspondentes ao invés de realizar aproximações locais em superfícies definidas localmente por planos tangentes, como proposto por Chen e Medioni (1992). Esta abordagem somente é possível, uma vez que após aplicar as etapas anteriores, restam nas nuvens de pontos LiDAR apenas pontos pertencentes a superfícies planas. Deste modo, não é necessário realizar aproximações locais, tendo em vista que há correspondência entre pontos e planos previamente extraídos.

O algoritmo proposto neste trabalho, também difere do método ICP original na etapa de seleção de pontos correspondentes. Para cada ponto da faixa LiDAR de referência (p1) é calculado um ponto pertencente ao plano da faixa LiDAR de pesquisa (p2). As coordenadas calculadas correspondem ao ponto do plano da faixa de pesquisa que é mais próximo do ponto da faixa de referência. A localização deste ponto no espaço é definida pela interseção do plano da faixa de pesquisa com a linha reta perpendicular ao plano que passa pelo ponto da faixa de referência, que representa o ponto correspondente à projeção do mesmo no plano da faixa adjacente. Basicamente, as coordenadas do ponto p2 são determinadas através da solução do seguinte sistema linear:

sendo p1 o ponto da faixa de referência, o vetor normal do plano da faixa de pesquisa (unitário), d2 a distância do plano da faixa de pesquisa à origem, p2 a projeção do ponto da faixa de referência no plano da faixa de pesquisa e s um escalar.

A primeira equação do sistema consiste na equação do plano da faixa LiDAR de pesquisa. A segunda equação corresponde à equação da linha reta perpendicular a esse plano que passa pelo ponto p1. Como p2 é a projeção de p1 no plano, a distância entre esses pontos corresponde à distância entre o ponto p1 e o plano definido pelos parâmetros e d2. Deste modo, ao aplicar na função erro do ICP os pontos p1 e a sua projeção p2, essa função passa a representar a minimização do quadrado da distância entre um ponto e seu respectivo plano. Com isso a referida equação toma a seguinte forma, a saber:

sendo o ponto i do conjunto b pertencente ao plano k e o ponto do conjunto a pertencente ao plano k que é mais próximo de .

Como há um numero significativamente inferior de planos em relação a pontos, os parâmetros de transformação são aplicados nos planos, resultando em melhor desempenho computacional. Após o cálculo dos parâmetros é possível avaliar os resíduos da distância ponto-a-plano antes e após a aplicação dos parâmetros calculados. Entretanto, não é possível avaliar as componentes altimétricas e planimétricas de modo independente. Esse processo é realizado na etapa seguinte empregando as linhas extraídas.

3.3.2 Avaliação Planimétrica e Altimétrica

As linhas de cumeada permitem determinar separadamente a precisão planimétrica e altimétrica dos dados LiDAR. Para uma linha reta que apresenta um ângulo de apenas 0,5º em relação à horizontal, por exemplo, o valor de h varia 1 cm a cada 11,45 m. Deste modo, pode ser considerado que uma linha de cumeada apresenta valores de h constantes ao longo de sua extensão.

A precisão altimétrica é obtida através da REMQ calculada a partir das altitudes de linhas de cumeada correspondentes, obtidas a partir da coordenada Z do centro aproximado da reta, como descrito na subseção 3.2. Já o cálculo da precisão planimétrica é baseado no trabalho de Vosselman (2008). Cada linha de cumeada extraída fornece a informação da componente planimétrica do deslocamento perpendicular à linha reta. Todos os cálculos são feitos no plano, sendo desprezada a coordenada Z. O cálculo do deslocamento planimétrico pode ser feito usando a Equação 6, como segue:

sendo e o erro planimétrico, xa e ya as coordenadas planimétricas do ponto pA1 de uma faixa, α2 e d2 os parâmetros da reta definida pela linha de cumeada correspondente (calculados a partir de p02 e ).

O erro planimétrico corresponde à distância entre o centro aproximado de uma linha e sua correspondente na faixa adjacente. Na ausência de erros, o ponto pA1 de uma faixa deveria pertencer à reta definida pela linha de cumeada extraída na faixa adjacente. A Figura 5 ilustra a obtenção do erro planimétrico.

Figura 5: Ilustração do erro planimétrico. 

A precisão planimétrica é obtida através da REMQ calculada a partir dos erros estimados pela Equação 6. A seguir serão apresentados os experimentos e resultados obtidos com o método proposto.

4. Experimento e análise dos resultados

Para avaliação do método proposto foi conduzido um experimento usando dados LiDAR referentes à varredura a LASER do bairro Bacacheri (Curitiba-PR) realizada em dezembro de 2006. A região é urbanizada e a área de sobreposição entre as duas faixas empregadas apresenta aproximadamente 2.000.000 m² (4 km x 500 m) sendo que cada uma apresenta densidade média de 1 ponto/m². O algoritmo desenvolvido para realização dos experimentos foi dividido em duas partes: a extração de feições foi implementada em C++ com auxílio da biblioteca Point Cloud Library; já o modelo matemático foi implementado na linguagem de programação Matlab.

4.1 Experimentos

O método proposto é automático, no entanto é necessário definir alguns parâmetros para o funcionamento eficiente e robusto dos algoritmos de filtragem, extração de linhas retas e correspondência entre as feições. Esses parâmetros foram empiricamente definidos para uma pequena parcela da área de estudo e, posteriormente, aplicados com êxito na totalidade dos dados. A Tabela 1 mostra os parâmetros definidos para cada uma das etapas do processo de filtragem, extração de linhas retas e correspondência entre as feições.

Tabela 1: Parâmetros definidos para cada uma das etapas do processo de filtragem, extração de linhas retas e correspondência entre as feições 

A Tabela 2 mostra a quantidade de pontos usados em cada etapa do método proposto.

Tabela 2: Quantidade de pontos em cada etapa. 

A Figura 6 apresenta, visualmente, os resultados obtidos após o emprego de das etapas de filtragem supracitadas.

Figura 6: Evolução do processo de filtragem dos dados: (a) Porção original da nuvem de pontos LiDAR; (b) Resultado após o emprego do filtro morfológico; (c) Resultado obtido após o emprego do filtro RANSAC; (d) Planos segmentados com o algoritmo de crescimento de regiões. 

A avaliação do método proposto foi dividida em duas partes. Na primeira parte, os resultados obtidos foram comparados considerando as seguintes situações:

- Situação ideal: os pontos da faixa de referência são aplicados sobre os planos extraídos da própria faixa, ou seja, esses valores seriam obtidos na hipótese dos pontos da faixa LiDAR de referência se adequarem perfeitamente aos planos da faixa LiDAR adjacente;

- Antes da transformação: os pontos da faixa LiDAR de referência são aplicados diretamente sobre os planos extraídos da faixa LiDAR adjacente sem qualquer transformação; e

- Método de Referência: solução obtida por ajustamento paramétrico não linear. Esta solução é baseada no método apresentado por Sande et al. (2010). Porém, é empregada a matriz de rotação associada a uma translação em vez de uma transformação afim 3D.

Já na segunda parte, foram empregados planos de verificação. Neste caso, foram aleatoriamente selecionados 25% dos planos extraídos. Por exemplo, na faixa de referência (Fx1**) foram extraídos 443 planos enquanto que na faixa de pesquisa (Fx2**) foram extraídos 560 planos. Entre os planos, foram encontrados 304 correspondentes. Também foram inicialmente extraídas 63 linhas de cumeeiras, das quais 15 foram descartadas pelo processo de eliminação. A Tabela 3 apresenta a distribuição de frequência da precisão dos planos com correspondências estabelecidas.

Tabela 3: Distribuição de frequência da precisão dos planos extraídos. 

A Tabela 4 apresenta a distribuição de frequência do número de pontos desses mesmos planos.

Tabela 4: Distribuição de frequência do número de pontos dos planos extraídos. 

Vale notar que os parâmetros para extração dos planos utilizados foram definidos buscando estabelecer uma relação ótima de quantidade versus qualidade de planos. A Tabela 5 apresenta os resultados obtidos após a realização da primeira parte da avaliação.

Tabela 5: Comparação dos resultados obtidos na 1ª parte da avaliação 

A Tabela 6 apresenta os parâmetros de translação e rotação (quatérnios convertidos para ângulos de Euler) calculados com o emprego do método proposto e pelo método de referência.

Tabela 6: Parâmetros obtidos na 1ª parte da avaliação 

A Figura 7a mostra o histograma da distância ponto-a-plano da situação ideal. Já na figura 7b pode ser visto, simultaneamente, os histogramas da distância ponto-a-plano antes (pior situação, na cor vermelha) e após a aplicação do método proposto (na cor azul). Os histogramas apresentam a quantidade de pontos no eixo vertical e a distância ponto-a-plano no eixo horizontal (em metros).

Figura 7 - Histogramas das distâncias ponto-a-plano: (a) situação ideal; (b) antes e após a aplicação dos parâmetros de transformação.

Figura 7: Histogramas das distâncias ponto-a-plano: (a) situação ideal; (b) antes e após a aplicação dos parâmetros de transformação. 

Após a aplicação da transformação, a avaliação da acurácia planimétrica e altimétrica é realizada através do cálculo da REMQ utilizando as 48 linhas retas correspondentes extraídas. O REMQ altimétrico encontrado é igual a 2,86 cm, enquanto que o REMQ planimétrico é igual a 12,60 cm. A Figura 8 mostra os histogramas dos erros altimétricos e planimétricos.

Figura 8: (a) Erros altimétricos; (b) Erros planimétricos. 

Em seguida, foi realizada a segunda parte da avaliação do método. Dos 304 planos, apenas 228 foram empregados no cálculo dos parâmetros de transformação, enquanto que 76 foram empregados no processo de verificação. A Tabela 7 apresenta os parâmetros de translação e rotação calculados pelo método proposto.

Tabela 7: Parâmetros calculados na 2ª avaliação 

Os resultados obtidos para a 2ª parte da avaliação do método proposto são apresentados na Tabela 8, onde foram avaliadas as distâncias ponto-a-plano antes e após a aplicação do método nos planos utilizados na verificação.

Tabela 8: Comparação dos resultados obtidos na 2ª avaliação 

Na Tabela 9 é apresentada uma comparação dos resultados obtidos na 2ª avaliação e na Figura 9 são apresentados a distribuição de frequências e o histograma da REMQ calculada por plano de verificação.

Tabela 9: Comparação dos resultados obtidos na 2ª parte da avaliação 

Figura 9: Histogramas da REMQ por plano de verificação antes e após a aplicação da transformação. 

A seguir serão apresentadas as discussões dos resultados obtidos.

4.2 Discussão dos Resultados

Como pode ser visto na Tabela 2 cada etapa do método proposto reduz significativamente o número de pontos de cada faixa analisada. A Tabela 3 mostra que o valor do desvio-padrão da maioria dos planos é em torno de 3 a 4 cm, sendo que quase a totalidade dos planos apresenta desvio-padrão entre 2 a 5 cm. A alta qualidade dos planos é consequência da baixa tolerância empregada na filtragem RANSAC e no refinamento realizado durante a extração de planos, que é de apenas 10 cm. O aumento da tolerância e a redução do número mínimo de pontos permitem o aumento do número de planos extraídos, mas em contrapartida aumentam a probabilidade da extração de falsos positivos ou também de planos imprecisos. Já na Tabela 4 pode ser verificado que 75% dos planos apresentam de 40 (tamanho mínimo) a 100 pontos. Na Tabela 5, os valores da média, desvio-padrão, REMQ e maior distância foram obtidos a partir das distâncias dos pontos em relação aos seus respectivos planos em cada um dos casos testados. Os valores de maior discrepância correspondem ao valor absoluto da maior distância entre ponto e plano encontrada em cada caso.

Analisando a coluna da situação ideal (ver Tabela 5) é possível verificar que os parâmetros dos planos foram corretamente determinados, tendo em vista que a média da distância entre os pontos e seus respectivos planos é um valor próximo de zero (restando apenas erros numéricos) e o seu histograma (ver Figura 7a) apresenta a forma de uma curva gaussiana. Também pode ser verificado o corte nos extremos (± 0,1 m) decorrente da tolerância aplicada durante o processo da extração de planos. Não foram encontradas tendências significativas, dado que o valor da REMQ é igual ao desvio padrão calculado. Em relação aos resultados obtidos com a pior situação, verifica-se a presença de erros sistemáticos, uma vez que a média da distância ponto a plano é aproximadamente 6 cm. Também pode ser verificado visualmente que o formato do histograma da distância ponto a plano se afasta de uma curva normal (em vermelho na figura 7b). Ainda na Tabela 5, os resultados obtidos com o método proposto mostram que a média das distâncias ponto-a-plano é -0,0019 cm, enquanto que a REMQ e o desvio-padrão são iguais a 6,04 cm e a maior discrepância encontrada foi de 26,10 cm. Além disso, o histograma gerado (em azul na figura 7b) apresenta a forma de uma curva gaussiana.

Comparando os resultados das Tabelas 5 e 6, verifica-se que os resultados do método de proposto são praticamente iguais ao do método de referência, tendo em vista que ambos minimizam a mesma função erro só que através de ferramentas matemáticas distintas. A pouca diferença pode resultar da implementação da aplicação das rotações, que no método proposto é realizada por quatérnios sendo que no método de referência são utilizados os ângulos de Euler.

Na Tabela 7, pode ser verificado que há uma translação significativa entre as faixas, uma vez que apenas a componente em X da translação é igual a 42 cm, indicando a presença de erros sistemáticos. Os valores encontrados para os parâmetros de translação indicam que as faixas empregadas no experimento não foram previamente ajustadas. Os valores dos ângulos de rotação calculados não foram significativos, apresentando valores inferiores a 0,02º.

Observa-se na Tabela 8 a redução significativa da REMQ nos dados empregados para verificação. Essa redução pode ser constatada também nos histogramas da Figura 10. Após a aplicação dos parâmetros calculados, o gráfico do histograma das distâncias ponto-a-plano assume a forma de uma curva gaussiana centrada em zero.

Figura 10: Histogramas de distância ponto-a-plano antes e após a aplicação da transformação nos dados de verificação. 

Como pode ser verificado na Tabela 9 e na Figura 10, antes da transformação aproximadamente metade dos planos apresentaram REMQ superior a 10 cm. Após a transformação, mais de 90% dos planos apresentaram REMQ inferior a 10 cm. A seguir serão apresentadas as conclusões e recomendações para trabalhos futuros.

5. Conclusões e recomendações

Este trabalho apresentou um método automático para avaliação da acurácia relativa de dados LiDAR aerotransportado. O método emprega simultaneamente primitivas extraídas dos telhados das edificações presentes nas faixas LiDAR obtidas por varredura LASER. Para investigar a potencialidade do método proposto, foi conduzido um experimento usando dados LiDAR de uma área urbanizada. Basicamente, três aspectos principais devem ser discutidos neste trabalho, isto é: o grau de automação do processo, a implementação de uma variação do ICP e a avaliação da acurácia relativa planialtimétrica dos dados LiDAR.

A combinação dos algoritmos RANSAC, Outlier Removal, filtro morfológico progressivo aliados ao algoritmo de crescimento de regiões resultou na extração automática de superfícies planas com precisão. A autonomia do método proposto propicia a redução do tempo de execução do projeto reduzindo, principalmente, a necessidade de intervenção humana. Além disso, o método independe de dados externos, devido ao fato de estar sendo estimada apenas a acurácia relativa, implicando na redução de custo na coleta de dados. Por outro lado, o pré-processamento e a filtragem dos dados agregaram um elevado custo computacional ao método.

As etapas que antecedem a determinação dos parâmetros de transformação corrigiram a deficiência apontada por Habib et al. (2010), que afirmaram que o ICP não apresenta um desempenho adequado na avaliação da qualidade devido a ausência de correspondências exatas, dado que a partir do método proposto para a etapa de filtragem é possível atribuir correspondências exatas entre pontos e planos. Deste modo, é possível minimizar as distâncias entre cada ponto da faixa de referência em relação ao seu respectivo plano da faixa de pesquisa.

A principal característica do modelo de otimização proposto é a possibilidade de avaliar detalhadamente a acurácia relativa entre faixas adjacentes, através do cálculo dos erros sistemáticos e análise das discrepâncias ponto-a-plano. Além disso, o emprego das linhas retas derivadas das cumeadas das edificações permitiu verificar simultaneamente a acurácia relativa planimétrica da altimétrica, sendo a grande contribuição do método proposto em relação ao estado da arte. Ainda é importante notar que o método proposto também pode ser empregado na avaliação da acurácia absoluta, desde que seja feita a correspondência entre os pontos de controle e superfícies planas encontradas nas nuvens de pontos 3D.

Como trabalhos futuros serão conduzidas investigações sobre a possibilidade do emprego do método proposto para avaliação da acurácia absoluta dos dados LiDAR, otimização computacional da etapa de filtragem dos dados e aperfeiçoamento do algoritmo de extração de linhas.

AGRADECIMENTOS

Os autores externam seus agradecimentos ao Exército Brasileiro pela concessão de licença remunerada para realização do curso de Mestrado em Ciências Geodésicas e ao CNPq por concessão de bolsa de Produtividade (no. Processo 302644/2013-0) ao orientador deste trabalho.

REFERÊNCIAS BIBLIOGRÁFICAS

López, Francisco Javier Ariza. Fundamentos de Evaluación de la Calidad de la Información Geográfica. 1. ed. Jaen: Universidade de Jaén, 2013. [ Links ]

Besl, Paul J.; Mckay, Neil D. A method of registration of 3-D shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence, v. 14, n. 2, p. 239-256, 1992. [ Links ]

Chen, Yung; Medioni, Gérard. Object modelling by registration of multiple range imagesImage and vision computing, v. 10, n. 3, p. 145-155, 1992. [ Links ]

Crombaghs, Marc; Brügelmann, Regine; de Min, Erik J. On the adjustment of overlapping strips of laser altimeter height data. International Archives of Photogrammetry and Remote Sensing, v. 33, n. B3/1, p. 230-237, 2000. [ Links ]

Csanyi Nora; Toth, Charles K.; Improvement of Lidar Data Accuracy Using Lidar-Specific Ground Targets. Photogrammetric Engineering & Remote Sensing, v. 73, n. 4, p.385-396, 2007. [ Links ]

Fischler, Martin A.; Bolles, Robert C. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, v. 24, n. 6, p. 381-395, 1981. [ Links ]

Friess, Peter. Toward a rigorous methodology for airborne laser mapping. EuroCOW, p. 25-27, 2006. [ Links ]

Habib, Ayman; Kersting, Ana Paula; Bang, Ki In; Lee Dong-Cheon. Alternative methodologies for the internal quality control of parallel LiDAR strips. IEEE Transactions on Geoscience and Remote Sensing, v. 48, n. 1, p. 221-236, 2010. [ Links ]

Horn, Berthold K. P. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America A, v. 4, n. 4, p. 629-642, 1987. [ Links ]

Rabbani, Tahir; Van Den Heuvel, Frank; Vosselmann, George. Segmentation of point clouds using smoothness constraint. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, v. 36, n.5, p. 248-253, 2006. [ Links ]

Sande, Corné Van Der; Soudarissanane, Sylvie; Khoshelham, Kourosh. Assessment of relative accuracy of AHN-2 laser scanning data using planar features. Sensors, v. 10, n. 9, p. 8198-8214, 2010. [ Links ]

Kilian, Johannes; Haala, Norbert; Englich, Markus. Capture and evaluation of airborne laser scanner data. International Archives of Photogrammetry and Remote Sensing, v. 31, p. 383-388, 1996. [ Links ]

Monico, João Francisco Galera; Dal Poz, Aluir Porfírio; Galo, Maurício; Dos Santos, Marcelo Carvalho; Oliveira, Leonardo Castro de. Acurácia e precisão: revendo os conceitos de forma acurada. Boletim de Ciências Geodésicas, v.15, n.3, 2009. [ Links ]

Pfeifer, Norbert C.; Elberink, Sander Oude; Filin, Sagi. Automatic tie elements detection for laser scanner strip adjustment., International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences v. 36, n. 3/ W3,p. 1682-1750, 2005. [ Links ]

Rentsch, Matthias; Krzystek, Peter. LiDAR strip adjustment using automatically reconstructed roof shapes. International Archives of the Photogrammetry, Remote Sensing and Spatial Information, 2009. p. 158-164, [ Links ]

Rusu, Radu Bogdan; Blodow, Nico; Beetz, Michael. Fast Point Feature Histograms (FPFH) for 3D Registration, IEEE International Conference on Robotics and Automation (ICRA), Kobe, Japan, May 12-17, 2009. [ Links ]

Skaloud, Jan; Lichti, Derek. Rigorous approach to bore-sight self-calibration in airborne laser scanning. ISPRS Journal of Photogrammetry and Remote Sensing, v. 61, n. 1, p. 47-59, 2006. [ Links ]

Toth, Charles K.; Paska, Eva; Brzezinska, Dorota. Using Road Pavement Markings as Ground Control for LiDAR Data., International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences v. 37, part B1, p.189-195, 2008. [ Links ]

Vosselman, George. Analysis of planimetric accuracy of airborne laser scanning surveys., International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences v. 37, part B3a, p. 99-104, 2008. [ Links ]

Vosselman, George. Automated planimetric quality control in high accuracy airborne laser scanning surveys., ISPRS Journal of Photogrammetry and Remote Sensing v. 74, p. 90-100, 2012. [ Links ]

Wehr, Aloysuys; Lohr, Uwe. Airborne Laser Scanning - An Introduction and Overview. ISPRS Journal of Photogrammetry & Remote Sensing. v. 54, pp. 68-82. 1999. [ Links ]

Zhang, Keqi; Chen, Shu-Ching; Whitman, Dean; Shyu, Mei-Ling; Yan, Jianhua; Zhang, Chengcui. A progressive morphological filter for removing nonground measurements from airborne LIDAR data. Geoscience and Remote Sensing, IEEE Transactions on, v. 41, n. 4, p. 872-882, 2003 [ Links ]

Received: September 2014; Accepted: August 2015

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License