Mestre em Sensoriamento Remoto

Universidade Federal de Minas Gerais

Instituto de Geociências - Departamento de Cartografia

Av Antônio Carlos, 6627 - Campus da Pampulha

Belo Horizonte - MG - CEP 31270-901 Fone (031) 4995428

Nas últimas décadas, o uso das imagens de satélites de sensoriamento remoto tem sido objeto de inúmeros estudos e pesquisas, visando a aplicações nas atividades de mapeamento e Sistemas de Informações Geográficas. As imagens de Sensoriamento remoto são obtidas em formato raster com dimensões limitadas que variam de acordo com o sistema sensor utilizado (MSS, TM, SPOT, AVHRR, etc). Quando se deseja analisar áreas extensas ou localizadas na vizinhança de duas ou mais imagens, surge a necessidade de juntá-las através de um mosaico.

Na montagem de um mosaico digital, a partir de duas ou mais cenas com faixas de superposição, a junção entre as imagens é geralmente aparente. As bordas aparentes são resultantes das mudanças ocorridas entre as épocas de aquisição das imagens. Neste intervalo de tempo, as condições atmosféricas podem ter variado, a umidade do solo pode estar diferente, mudanças sazonais podem ter ocorrido na vegetação, problemas de calibração do sensor podem ter acontecido, etc. Freqüentemente, as bordas criadas artificialmente na região de junção são mais visíveis que a informação verdadeira contida na imagem, e podem provocar erros de interpretação.

Assim, o objetivo do presente trabalho é apresentar uma métodologia para construir mosaicos de imagens digitais de sensoriamento remoto que solucione o problema apresentado. A metodologia utiliza uma abordagem que contempla a transformação das imagens para um sistema de projeção cartográfica, modela a transformação geométrica através de polinômios com uso de pontos de controle de terreno, reamostra as imagens, ajusta os seus níveis de cinza para um padrão comum, escolhe a linha irregular mais apropriada para junção das imagens e suaviza a transição entre elas nas áreas próximas à linha de junção utilizando um filtro linear de efeito degradê

No enfoque deste trabalho, a montagem de mosaicos está fundamentada nas linhas gerais definidas de forma simplificada por Milgram (1975) e refinadas em Milgram (1977). Contribuiram também os trabalhos desenvolvidos por Zobrist et al (1983), Hood et al (1989), Miller (1989) e Shiren et al (1989). Partindo destas linhas bases foram acrescentados no trabalho vários refinamentos, dentre os quais o detalhamento e algumas variações adaptativas dos algoritmos, a implementação em liguagem C, e a testagem da efetividade dos algoritmos em imagens do sensor Thematic Mapper (TM) do satelite LANDSAT visando aplicações nos trabalhos de mapeamento e Sistemas de Informações Geográficas. Em resumo, quatro estágios são necessários na montagem de um mosaico digital sem junções aparentes, segundo a abordagem deste trabalho.

1) O primeiro passo requer o perfeito registro das imagens. Geralmente são utilizando dois processos: o modelo inverso da geometria de imageamento ou a modelagem por polinômios que, neste caso, foi a escolhida e será discutida mais adiante.

2) No segundo passo, é necessário o ajuste da tonalidade de cinza de uma imagem em relação a outra, de forma que a média dos níveis de cinza nas regiões de superposição sejam semelhantes para ambas. Para este proposito geralmente utilizam-se técnicas de ajuste de histograma (Richards, 1986). Após este passo as bordas artificiais tornam-se significativamente menos visíveis, apesar de ainda serem aparentes.

3) No terceiro passo, é criada uma linha irregular para junção das imagens. Um ponto de junção é escolhido para cada linha. A informação de cada linha no mosaico final à esquerda do ponto de junção será puxada do segmento de linha da imagem esquerda, e a imagem direita fornecerá a informação para o restante da linha do mosaico. O ponto de junção será escolhido onde o somatório das diferenças de níveis de cinza sobre uma região de w pontos é mínimo, isto é, onde o valor da função Df(n) abaixo é mínimo.

onde,

E(n) (n=1,2,...,K) são os valores dos pixels da linha considerada da imagem esquerda na região de superposição;

D(n) (n=1,2,...,K) são os valores dos pixels da linha considerada da imagem direita na região de superposição; n é um índice que indica a posição do pixel pesquisado;

K é a largura da região de pesquisa que terá que ser menor do que a região de superposição; e w é uma máscara que percorre toda a extensão K, correspondente à largura da região para cálculo do somatório das diferenças (Df(n)), que deve ser menor ou igual a K. No intuito de evitar o surgimento de bordas horizontais, um ótimo refinamento para este passo consiste em restringir o intervalo de escolha do ponto de junção de uma dada linha, dependendo da magnitude do Df(n) mínimo da linha anterior. Se o ponto de junção anterior teve um Df(n) mínimo elevado, é conveniente escolher o ponto de junção atual próximo do anterior. Se o Df(n) mínimo foi pequeno, a posição do ponto de junção não será tão crítica. Este refinamento incorpora a grande vantagem de contornar nuvens e outras feições presentes numa imagem e ausentes na outra.

4) Finalmente, o quarto passo consiste na suavização de eventuais abruptas diferenças de níveis de cinza na vizinhança do ponto de junção que por qualquer motivo tenham permanecido. Isso pode ser obtido através de uma transição gradual ou efeito degradê.

Precedendo a montagem do mosaico digital algumas etapas preparatórias tornaram-se indispensáveis, uma vez que o objetivo final é o seu uso em trabalhos de mapeamento e Sistemas de Informações Geográficas devendo, portanto, obedecer a critérios e padrões de precisão já consagrados pelo uso e até mesmo estabelecidos na legislação cartográfica própria de cada país.

Vários estudos cientificos já demonstraram que as imagens do sensor TM com correções geométricas padrões de sistema possuem excelente precisão interna (Bryant et al 1985, Borgerson et al 1985, d'Alge 1987). Entretanto, para efeito de mapeamento e aplicações em SIG as imagens precisam ser registradas com Pontos de Controle de Terreno (PCT) por três razões básicas a saber.

1) Eliminar os erros grosseiros de posicionamento, ou seja os erros de translação em latitude e longitude apresentados pelas imagens com correção de sistema sem uso de PCT (baseada somente em dados de efemérides).

2) Orientar os eixos das imagens com relação ao Norte e Este geográficos, tendo em vista que seus eixos estão originalmente orientados com relação à linha de vôo do satélite e a direção da varredura do sensor.

3) Estabelecer o sistema de projeção cartográfica desejado, neste caso o Sistema Universal Transverso de Mercator (UTM), e fornecer uma referência comum para permitir o mosaico das imagens.

Para cumprir este objetivo é necessário que sejam definidos Pontos de Controle que tenham coordenadas no sistema de cada imagem e no sistema de projeção UTM, calcular coeficientes de transformação entre os sistemas de cada imagem e o sistema UTM e reamostrar cada imagem no sistema UTM.

A forma mais confiável de obtenção de pontos de controle é através de levantamento de campo, utilizando métodos geodésicos ou rastreamento de satélites de posicionamento. Entretanto, este é um processo custoso em termos de tempo e dinheiro, estando ao alcance de poucas instituições e empresas ou projetos de maior porte. Como alternativa menos adequada, porém mais barata, utilizam-se mapas em escalas grandes como fonte de pontos de controle. No caso deste trabalho, por limitações de custos, utilizou-se este processo. Os PCT foram extraídos das folhas de Mangaratiba e Angra dos Reis, na escala de 1:50.000 do Mapeamento Sistemático Nacional.

Neste estudo, foram selecionados vários pontos de controle para cálculo dos coeficientes de transformação entre os sistemas imagem-mapa, e outra quantidade similar para avaliação da qualidade da transformação, através da analise do erro médio quadrático dos resíduos (EMQ), conforme pode ser visto nas Tabelas 1 a 4, apresentadas na seção ANÁLISE DOS RESULTADOS. A expressão do EMQ é dada por

Onde, X e Y são as posições no sistema de projeção do mapa ou terreno conhecidas a priori X" e Y" são as posições correspondentes, calculadas utilizando os coeficientes de transformação resultantes do modelo e n é o numero de PCT utilizados na avaliação.

Durante o processo de seleção é inevitável a rejeição, passo a passo, dos pontos que apresentarem resíduos grosseiros, até que se atinja um EMQ compatível com o Padrão de Exatidão Cartografica (PEC) dos mapas nacionais. O PEC para os documentos cartográficos classe A estabelece que 90% dos pontos contidos no mapa devem estar dentro de uma incerteza de 0.5 milímetros na escala do mapa.

Considerando imagens TM processadas com correção de sistema que apresentam boa precisão interna, escolheu-se um modelo polinomial de primeiro grau para a transformação geométrica das imagens (Equações 2), conforme recomenda Welch et al. (1985). Este modelo corrige uma rotação, uma translação no eixo X, outra translação no eixo Y e um fator de escala para cada um destes eixos, mantendo, porém, o paralelismo entre eles. As equações do modelo são dadas por

colocando-se na forma matricial tem-se

L = A.B          (3)

cuja solução para os coeficientes é dada por (Neter e Wasserman, 1974)

B = (A^TA)^-1 . A^TL         (4)

onde B é o vetor de coeficientes de transformação, A é a matriz com as coordenadas de imagem e L é o vetor de coordenadas de terreno no sistema UTM

Com base na equação 4 foi desenvolvido o programa CALCULA-COEFICIENTES, em linguagem C que recebe como entrada um arquivo de Pontos de Controle e fornece um arquivo de saída contendo quatro conjuntos de seis coeficientes de transformação. Para um detalhamento mais completo sobre os algoritmos e a codificação dos programas em liguagem C consultar Elmiro (1994). Determinados os coeficientes de transformação, o passo seguinte é a reamostragem das imagens no sistema UTM. Foi desenvolvido um módulo (REAMOSTRAGEM) que lê o arquivo de coeficientes, as imagens esquerda e direita, reamostra as duas imagens utilizando um interpolador de vizinhança mais próxima e fornece como saída: as imagens esquerda e direita corrigidas e acompanhadas de um arquivo contendo dados de registro.

O módulo MOSAICO foi programado para identificar a área de superposição das imagens, calcular a média de brilho de ambas nesta área e ajustar a tonalidade da imagem da direita para a média da imagem esquerda. Após o ajuste das médias, o módulo cria automaticamente uma linha irregular para junção das imagens, baseado na menor diferença radiométrica local (Equação 1). O módulo aplica também um filtro de efeito degradê para transição das imagens nas proximidades do ponto de junção. O filtro elimina eventuais diferenças bruscas de tonalidade que não tenham sido corrigidas pelo ajuste de médias. Assim, cada ponto Pi de uma da vizinhança de v pixels (sendo (v-1)/2 pixels de cada lado do ponto de junção) será calculado por

onde,

i é um índice relativo que varia de 1 a v, sendo o ponto de junção, o ponto central

Pi é o nível de cinza calculado para o pixel de índice i da vizinhança.

Ei é o nível de cinza do pixel da imagem esquerda de índice i da vizinhança

Di é o nível de cinza do pixel da imagem direita de índice i da vizinhança

Obtém-se, assim, para cada ponto uma média onde os pesos são ponderados, de acordo com o lado e a distância do ponto de junção. Em pseudo-linguagem de programação o módulo MOSAICO terá a lógica abaixo, a codificação em linguagem C pode ser encontrada em Elmiro (1994) .

Chamar módulo para ajuste_radiométrico

Começar pela imagem mais alta

ENQUANTO pixel processado não atingir a imagem mais baixa

Cada linha do mosaico recebe zeros no lado referente à imagem mais baixa mais a linha da imagem mais alta

ENQUANTO pixel processado não atingir o final de uma das imagens

ENQUANTO pixel da linha não atingir área de superposição

Linha do mosaico recebe linha da imagem esquerda

Chamar módulo para cálculo_do_ponto_de_junção

ENQUANTO pixel da linha dentro da área de superposição

SE linha não atingiu ponto de junção

Linha do mosaico recebe linha da imagem esquerda

SE linha passou do ponto de junção

Linha do mosaico recebe linha da imagem direita

ENQUANTO pixel da linha não atingir o final da linha

Linha do mosaico recebe linha da imagem direita

Chamar módulo suavização_da_junção

ENQUANTO pixel processado não atingir o final da última imagem

Cada linha do mosaico recebe zeros no lado referente à imagem que já terminou mais a linha da imagem que não acabou

O módulo para ajuste_radiométrico tem a forma abaixo

Começar pela imagem mais alta

ENQUANTO não atingir a imagem mais baixa

Prosegue leitura

ENQUANTO não atingir o fim de uma das imagens

ENQUANTO linha não atingir área de superposição

Prosegue leitura

ENQUANTO linha dentro da área de superposição

Acumular pixel da imagem esquerda

Acumular pixel da imagem direita

Calcular média da imagem esquerda

Calcular média da imagem direita

Calcular diferença de médias

Ajustar diferença de média a todos os pixels da imagem direita

A forma geral do módulo para cálculo_do_ponto_de_junção é

Receber largura da área de procura K e da máscara w

Começar no início da área de procura

ENQUANTO linha dentro da área de procura

Calcular somatório das diferenças Df(n) para os w pixels vizinhos

Armazenar cada valor Df(n) e a correspondente posição n dentro da linha

Calcular valor Df(n) mínimo

Retornar posição n do valor Df(n) mínimo para o módulo mosaico

O módulo para suavização_da_junção tem a forma

VARIANDO i de 1 a v em torno do ponto de junção

Nível de cinza do mosaico Pi = ((v-i) * nível_de_cinza-da _imagem_esquerda + i*nível_de_cinza_da_imagem_direita)/v

FIM-VAR

O enfoque dado neste trabalho à montagem de mosaicos limitou-se à utilização de duas imagens vizinhas. Para o caso de três ou mais imagens vizinhas é suficiente que se façam algumas adaptações aos algoritmos. No caso da junção das imagens, não haverá grandes dificuldades, sendo suficiente que o programa para cálculo dos coeficientes forneça também dados de registro relativos às demais imagens para que estes dados, por sua vez, sejam utilizados pelos demais programas para reamostragem e junção. No caso do ajuste de tonalidade, todas as imagens poderão ter suas médias inicialmente ajustadas a uma das imagens que servirá de referência radiométrica. Outra solução possível é mosaicar as imagens duas a duas de modo que o mosaico resultante das duas primeiras seja uma imagem de entrada para um mosaico com a terceira. Este resultado, por sua vez, será também uma imagem de entrada para um mosaico com a quarta imagem e assim sucessivamente.

A finalidade deste passo é verificar se o mosaico preserva a precisão geométrica das imagens individuais e atinge o Padrão de Exatidão Cartográfica (PEC) dos mapas nacionais classe A. Os testes devem ser feitos através de uma transformação de similaridade ou isogonal, utilizando um conjunto de Pontos de Controle para cálculo dos coeficientes e outro conjunto para avaliação do EMQ dos resíduos. Esta transformação não introduz deformações, sendo adequada para verificação do grau de similaridade do mosaico final com relação ao mapa e as imagens componentes (Silva et al., 1987).

A formulação matemática para a transformação isogonal é dada por

cuja solução para os coeficientes na forma matricial é dada por

B = (A^TA)^-1 . A^TL              (7)

Nas Tabelas 1 a 4 são apresentados os resíduos e o EMQ dos pontos de controle utilizados nas transformações geométricas executadas sobre as imagens para montagem do mosaico. Os dados foram obtidos através do programa CALCULA-COEFICIENTES. Pelos resultados observa-se que o EMQ, tanto para o cálculo dos coeficientes como para avaliação da qualidade da transformação é de ordem inferior a um pixel do sensor TM (30 metros). Estes resultados reforçam a teoria sobre a excelente precisão geométrica interna das imagens TM com correções geométricas de sistema.

TABELA 1 - RESÍDUOS DOS PCT UTILIZADOS PARA O CÁLCULO DOS COEFICIENTES NA IMAGEM ESQUERDA DO MOSAICO (valores em metros)

EMQ E = 10.16 EMQ N = 13.44

EMQ total = 16.85

TABELA 2 - RESÍDUOS DOS PCT UTILIZADOS PARA AVALIAÇÃO DA TRANSFORMAÇÃO POLINOMIAL NA IMAGEM ESQUERDA DO MOSAICO (valores em metros)

EMQ E = 7.42 EMQ N = 16.48

EMQ total = 18.07

TABELA 3 - RESÍDUOS DOS PCT UTILIZADOS PARA O CÁLCULO DOS COEFICIENTES NA IMAGEM DIREITA DO MOSAICO (valores em metros)

EMQ E = 8.34 EMQ N = 18.97

EMQ total = 20.72

TABELA 4 - RESÍDUOS DOS PCT UTILIZADOS PARA AVALIAÇÃO DA TRANSFORMAÇÃO POLINOMIAL NA IMAGEM DIREITA DO MOSAICO (valores em metros)

EMQ E = 12.69 EMQ N = 11.81

EMQ total = 17.34

A Figura 1 apresenta um pedaço de imagem da banda TM4, extraída do quadrante B, órbita ponto 218/76, obtida no dia 11/set/90 e a Figura 2 mostra um pedaço de imagem da banda TM4, quadrante A, órbita ponto 217/76 (orbita vizinha) obtida em 24/fev/90.

As Figuras 3 e 4 apresentam as imagens anteriores reamostradas para o sistema de projeção UTM.

A Figura 5 é um mosaico da banda TM4, com 946 colunas x 587 linhas, sem correção obtido a partir das duas imagens das Figuras 3 e 4, utilizando um programa de mosaico sem correção elaborado para testes do trabalho, consultar Elmiro (1994) para maior detalhamento. Observa-se que as diferenças radiométricas são significativas e a linha reta de junção entre as imagens é notória, tornando-se um resultado inadequado para fins práticos.

A Figura 6 apresenta um mosaico de dimensões semelhantes ao anterior a partir das mesmas imagens das Figuras 3 e 4. Neste caso porém, a média radiométrica da imagem direita foi ajustada para a média da imagem esquerda aplicou-se também uma linha de junção irregular para as duas imagens com base na menor diferença radiométrica local, Df(n), definida na Equação 1. Neste caso a largura da região de pesquisa foi estabelecida em K = 40 pixels e o tamanho da máscara de integração em w  = 10 pixels. Uma boa forma de se chegar a resultados satisfatórios é fazer testes interativos variando estes dois parâmetros. Utilizou-se aqui o programa MOSAICO com correção.

O mosaico final foi avaliado através da análise dos resíduos de uma transformação isogonal (Equação 6) e apresentou bons resultados, os quais podem ser visualizados nas Tabelas 5 e 6. O EMQ de 17 metros, resultante da avaliação mostra que o mosaico conserva uma qualidade geométrica compatível com as imagens componentes e se mantém dentro do PEC dos mapas nacionais classe A para escalas menores que 1:100.000. A escala de 1:100.000 permite um um EMQ de até 50 metros (0.5 mm x 100.000). Visualmente não se consegue identificar descontinuidades no mosaico, mesmo utilizando ampliações nas regiões costeiras e ilhas onde os detalhes são precisos e muito bem definidos.

A Figura 7 apresenta uma ampliação por interpolação linear de um trecho do mosaico na região da junção onde se pode ver com maior clareza a otima qualidade do trabalho. A Figura 8 mostra a composição colorida do mosaico das bandas TM5, TM4 e TM3 em vermelho, verde e azul respectivamente.

TABELA 5 - RESÍDUOS DOS PCT UTILIZADOS PARA O CÁLCULO DOS COEFICIENTES NA AVALIAÇÃO DO MOSAICO FINAL (valores em metros)

EMQ E = 12.59 EMQ N = 11.97

EMQ total = 17.37

TABELA 6 - RESÍDUOS DOS PCT UTILIZADOS PARA AVALIAÇÃO DO MOSAICO FINAL (valores em metros)

EMQ E = 12.34 EMQ N = 12.71

EMQ total = 17.71

Para testar a efetividade dos algoritmos a região em estudo foi propositalmente escolhida em condições adversas, ou seja com relevo bastante movimentado, muitas ocorrências de sombras topográficas e as imagens vizinhas tomadas em diferentes estações do ano. Fazendo uma análise geral da boa qualidade dos mosaicos obtidos e considerando as condições adversas do teste é de se esperar que resultados ainda melhores possam ser atingidos em regiões planas e utilizando imagens da mesma estação do ano, onde haverá menor ocorrência de sombras topográficas, menores mudanças sazonais, deslocamentos de relevo desprezíveis e onde a tonalidade radiométrica das imagens é mais uniforme e por cosequência as descontinuidades iniciais a serem corrigidas serão sensivelmente menores.

Fig. 1 - Imagem da banda TM4 quadrante B, componente esquerda do mosaico.

 

 

Fig. 2 - Imagem da banda TM4 quadrante A, componente direita do mosaico.

 

 

Fig. 3 - Imagem da Figura 1 reamostrada para o sistema UTM.

 

 

Fig. 4 - Imagem da Figura 2 reamostrada para o sistema UTM.

 

 

Fig. 5 - Mosaico da banda TM4, sem correção.

 

 

Fig. 6 - Mosaico da banda TM4 com correção apresentando ótimos resultados.

 

 

Fig. 7 - Detalhe ampliado do mosaico da Figura 6 na área de junção.

 

 

Fig. 8 - Composição colorida das bandas TM5, TM4 e TM3 em RGB.