Fundamentos Geoestatísticos

Referência principal

Aumond et al. (2018) — Kriging-based spatial interpolation from measurements for sound level mapping in urban areas. JASA 143(5), pp.2847-2857.

Relação densidade × desempenho (Aumond, Fig. 9)

Nosso caso: Belo Horizonte

• Área: 331 km²
• Medições: 7.450 no total, distribuídas em 10 cenários temporais (670–788 por cenário)
• Coleta: controlada e dedicada, com equipamento Classe 1 e dataloggers segundo a segundo
• Densidade efetiva por cenário: ~2.2 pts/km²
• RMSE obtido (RK Híbrido, trend-only): 5.76–8.02 dB
• RMSE obtido (IDW baseline): 6.64–8.96 dB

Referências complementares:

• Diretiva EU 2002/49/EC — tolerância de ±4.6 dB para modelos CNOSSOS-EU
• Can, Dekoninck & Botteldooren (2014) — RMSE 2–4 dB em escala de bairro
• Harman, Koseoglu & Yigit (2016) — comparação IDW/Kriging em Isparta, Turquia

Observado no projeto

• Nugget ratio (variograma global): 0.79 (nugget=36.3, sill=9.8, total=46.1 dB²)
• Nugget ratio RK resíduos: 0.83–0.91

Interpretação segundo a literatura

Aumond et al. obtiveram nugget = 0 com ~1500 pts/km² (referência densa). O nugget representa variabilidade não explicada pela estrutura espacial:

1. Variabilidade micro-escala — diferenças de ruído entre pontos muito próximos (fachadas, recuos, vegetação) que a amostragem esparsa não resolve
2. Ruído de medição — incerteza instrumental, eventos sonoros transitórios
3. Subamostragem — com 2.2 pts/km², o variograma empírico não consegue estimar a semivariância em curtas distâncias

Por que o nugget ratio dos resíduos RK é ainda maior (0.83–0.91)?

A regressão (trend) remove a variabilidade explicável pelas covariáveis. O que sobra (resíduos) é majoritariamente ruído aleatório sem estrutura espacial. Isso é esperado quando a tendência já captura a maior parte do sinal.

Implicação prática: Quando nugget_ratio > 0.8, o kriging dos resíduos não melhora a predição. O modelo opera em regime trend-only: a predição é inteiramente determinada pela regressão sobre covariáveis, sem componente geoestatística nos resíduos.

Regime trend-only: evidência empírica em BH

Nos 10 cenários do projeto, o compare_models (validação cruzada LOO interna do variograma) rejeitou todos os modelos de variograma em todos os cenários — nenhum modelo (exponencial, esférico, gaussiano) produziu erros padronizados com variância próxima de 1.0 e média próxima de 0. O fallback exponencial resulta consistentemente em nugget_ratio ≥ 0.80, ativando o regime trend-only.

Isso significa que, com a densidade amostral atual (~2.2 pts/km²), o Regression Kriging opera de fato como uma Regressão Espacial pura. O nome “Regression Kriging” descreve o framework metodológico (regressão + krigagem dos resíduos), mas na prática a krigagem dos resíduos contribui com zero — é tentar força-la degrada catastroficamente os resultados (cenário 8: RMSE explodiu de 6.0 para 25.8 dB antes da correção).

O ganho do RK sobre o IDW vem inteiramente da regressão com covariáveis:

Redução média: 12.7% no RMSE, com R² passando de ~0.01 (IDW) para ~0.25 (RK).

Observado vs. literatura

Diagnóstico

Nossos ranges elevados são provavelmente artefato de dois fatores:

1. Maxlag muito grande — o variograma empírico é calculado até o extent completo da cidade (~20 km), incluindo pares de pontos a distâncias onde a semivariância já estabilizou
2. Baixa densidade — poucos pares de pontos em curtas distâncias tornam os primeiros lags do variograma instáveis

Recomendação

Limitar maxlag do variograma empírico para ~2000m (≈ distância média entre medições). Isso deve revelar melhor a estrutura espacial de curto alcance e produzir ranges mais realistas (~500–2000m).

Justificativa

• Aumond et al. usam Matérn com κ=0.5, que é matematicamente equivalente ao modelo exponencial
• O modelo exponencial é o mais robusto para dados ambientais com transições suaves (Webster & Oliver, 2007)
• Nosso auto-fit testa exponential, spherical e gaussian — exponential é consistentemente selecionado como fallback

Evidência no artigo

• Universal Kriging (UK) com 4 distâncias a categorias de vias reduz o sill pela metade (22.4 → 11.2 dB²)
• A tendência melhora mais o resultado do que trocar distância euclidiana por distância na rede viária
• Conclusão dos autores: “efforts should be focused on the trend definition”

Nosso modelo RK Híbrido

Usa 14 covariáveis (distâncias viárias, transporte público, morfologia urbana, amenidades por perfil de atividade):

• Melhoria RK vs IDW: 10–15% redução no RMSE (média 12.7%)
• Sinal consistente: covariáveis selecionadas mantêm sinal estável entre 10 cenários

Categorias de vias (Aumond, Table 1 vs. nosso modelo)

Amenidades por perfil de atividade (abordagem original)

A covariável genérica osm_n_amenidades_250m (contagem total de amenidades num raio de 250m) foi substituída por duas covariáveis com perfil temporal distinto:

Justificativa: Os 10 cenários do projeto cobrem faixas horárias distintas em dois perfis de dia (ter/qua e sex/sáb). A contagem genérica não distingue áreas com concentração de bares (ruidosas à noite) de áreas com concentração de comércio (ruidosas de dia). A separação permite que o modelo RK atribua coeficientes diferentes por cenário.

Critérios de classificação:

• Postos de combustível (fuel) excluídos do comércio — o ruído associado é veicular, já capturado pelas covariáveis viárias
• osm_dist_hospital mantida como covariável separada de distância — hospitais são zonas de silêncio por legislação (Lei 9.505/08), não fontes de ruído

Achado de Aumond et al.

• Distância pela rede viária (ROAD) melhora OK, mas não melhora UK/RK
• Quando o trend já inclui distâncias a categorias de vias, usar ROAD é redundante

Implicação para o projeto

Nosso RK já incorpora distâncias viárias como covariáveis na regressão. Portanto, distância euclidiana é suficiente para o variograma dos resíduos — a informação viária já está no trend.

Abordagem do artigo

• Bootstrap com 1000 replicações de subconjuntos aleatórios
• Comparação contra mapa de referência denso (4360 pontos)

Nossa abordagem

• LOO Cross-Validation (Leave-One-Out): remove 1 ponto, prediz com os N-1 restantes, repete N vezes
• Métricas: RMSE, MAE, R², ME (bias), erros padronizados

Limitações conhecidas

• LOO-CV é pessimista com poucos dados (cada fold remove ~0.13% da informação)
• R² negativo em LOO-CV indica que a predição é pior que a média — esperado com nugget dominante em OK
• O bias próximo de zero (~0.0 dB) indica que o modelo não tem viés sistemático, mesmo com RMSE alto

Observado

• IDW (baseline): R² = -0.013 a 0.050 (varia por cenário)
• RK Híbrido (trend-only): R² = 0.19 a 0.29

O que significa

R² em LOO-CV compara o modelo contra a média global. R² < 0 significa que predizer a média seria melhor que o interpolador para aquele subconjunto.

Isso é esperado para IDW em escala municipal — sem covariáveis, o IDW com ponderação inversa converge para a média local, que pode ser pior que a média global em certas configurações espaciais.

RK melhora o R² consistentemente porque a tendência (regressão) captura variância que o IDW não consegue.

O modelo é adequado para o contexto?

Sim, com as seguintes ressalvas:

1. A limitação principal é a densidade amostral (~2.2 pts/km² por cenário), não a metodologia. O RMSE obtido é coerente com o esperado pela literatura para essa densidade. Note-se que o projeto totaliza 7.450 medições controladas com equipamento Classe 1, mas distribuídas em 10 cenários temporais.
2. O RK Híbrido é a melhor escolha disponível — a tendência com covariáveis urbanas é o componente que mais contribui para a qualidade. Isso está alinhado com a Conclusão de Aumond et al.
3. O viés próximo de zero (ME ≈ 0.0 dB em todos os cenários) indica que o modelo é não-enviesado, mesmo com RMSE moderado. Isso é a métrica mais importante para zoneamento.
4. A melhoria RK vs IDW (8–15% no RMSE) demonstra que as covariáveis capturam informação real, não ruído.
5. Evolução em relação ao estudo anterior (Martins & Murta, 2012): os mesmos autores avançaram de medições aproveitadas de ações fiscais (ambiente interno, ruído de fundo) para coleta controlada em campo aberto (ruído total L_eq), com metodologia de cenários e equipamento Classe 1.

O que não pode ser melhorado sem mais dados

• O RMSE não vai abaixo de ~5 dB com 2.2 pts/km² — isso é limite teórico
• O nugget ratio não vai diminuir significativamente — é reflexo da amostragem
• O R² não vai acima de ~0.30 em escala municipal — variância total muito alta

O que pode ser melhorado com o modelo atual

• Ajustar maxlag do variograma para ~2000m (pode melhorar marginalmente)
• Testar modelo Matérn com κ=0.5 explicitamente (equivalente ao exponencial, mas com ajuste mais flexível)
• Incorporar distância pela rede viária para o variograma OK (melhoria esperada modesta)

Descrição metodológica

O Monitoramento do Ambiente Sonoro Urbano de Belo Horizonte foi realizado com base em medições de campo, seguindo critérios técnicos e normas nacionais. As medições foram feitas em 3 (três) períodos ao longo de 3 (três) anos: de julho a novembro de 2023, março a abril de 2024 e setembro de 2024 a fevereiro de 2025. Em cada um dos pontos, a medição teve duração de 5 (cinco) minutos, tempo considerado suficiente para um diagnóstico preliminar da paisagem sonora da cidade.

Equipamento e normas

Os dados foram coletados com sonômetros da marca Svantek, modelo SV 971A, equipamentos de classe 1 conforme as normas internacionais, devidamente calibrados em laboratório pela última vez em março de 2023. As medições seguiram os critérios estabelecidos pela ABNT NBR 10.151 (versão atualizada de 2020), que dispõe sobre medição e avaliação de níveis de pressão sonora em áreas habitadas. Considerando a diversidade dos cenários urbanos — de vias alargadas a becos estreitos — buscou-se respeitar a distância mínima de 2 metros de paredes, muros, veículos ou outros objetos que pudessem refletir as ondas sonoras.

Faixas horárias e dias de medição

As medições foram realizadas em 5 faixas de horário:

Os dias de medição definidos foram terça-feira e quarta-feira (representando dias úteis típicos) e sexta-feira e sábado (representando dias de fim de semana). A escolha dessas combinações se deve ao fato de que, principalmente no período noturno, os padrões de tráfego nesses pares de dias são similares, favorecendo uma leitura mais consistente do ambiente sonoro.

Malha de pontos de medição

Os pontos de medição foram definidos por meio de uma malha quadricular sobre o mapa da cidade, com quadrados de aproximadamente 400×400 metros — e de 200×200 metros na região Centro-Sul, no perímetro dentro da Avenida do Contorno — posicionados nos vértices para cobrir uniformemente o território urbano.

De acordo com o cenário apresentado, tornou-se necessário adaptar a unidade-dia dos pontos de medição dessa malha, para uma análise alternada da convergência desses pontos, agrupadamente, distribuídos entre os dias típicos (terça e quarta-feira) e os de fim de semana (sexta e sábado).

Alguns pontos também precisaram ser realocados — por exemplo, quando coincidiam com pontos de ônibus, semáforos ou cruzamentos, locais onde o ruído é influenciado por frenagens, acelerações ou concentração de veículos. Nesses casos, os pontos foram deslocados apenas por alguns metros.

Contagem de veículos

A partir do terceiro período de monitoramento (setembro de 2024 a fevereiro de 2025), incluiu-se a contagem de veículos que passavam por cada local durante os 5 minutos de medição, com foco nos horários diurnos. Os veículos foram classificados em:

• Leves: carros
• Pesados: ônibus, caminhões, kombis, carros-fortes, vans
• Motocicletas

Essa classificação permitirá, em análises futuras, compreender melhor a influência do tráfego no ambiente acústico da cidade.

Visão geral

O modelo RK Híbrido utiliza 14 covariáveis espaciais para capturar a tendência (trend) do ruído urbano. Cada covariável associa a cada célula da grade (10×10m) um valor numérico que descreve uma característica do entorno — distância a vias, concentração de atividades, morfologia urbana etc. A regressão aprende qual combinação linear dessas variáveis melhor prevê o nível de ruído.

Todas as covariáveis são calculadas em Python (scripts/etapas/covariaveis.py) usando dois tipos de operação espacial:

1. Distância mínima — distância euclidiana da célula ao ponto mais próximo de um conjunto (ex: distância ao ponto de ônibus mais perto)
2. Contagem em raio — quantos pontos de um conjunto existem dentro de um raio fixo (ex: número de pontos de ônibus em 250m)

Ambas operações dependem de busca eficiente entre milhões de pares coordenada–coordenada. A estrutura de dados que viabiliza isso é a cKDTree.

O que é cKDTree

A cKDTree (C implementation of K-Dimensional Tree, da biblioteca scipy.spatial) é uma árvore binária de particionamento espacial otimizada para buscas de vizinhança em coordenadas multidimensionais.

Problema que resolve: dada a grade de BH com ~3.3 milhões de células e, por exemplo, ~11.000 pontos de ônibus, calcular a distância mínima de cada célula ao ponto de ônibus mais próximo exige, por força bruta, 3.3M × 11K = 36 bilhões de cálculos de distância. Na prática: inviável.

A cKDTree resolve isso em tempo O(n · log m) em vez de O(n · m), onde n é o número de células e m o número de pontos de referência. Para BH, isso reduz de ~36 bilhões para ~50 milhões de operações — uma redução de ~700×.

Como funciona

A árvore organiza os pontos de referência (ex: pontos de ônibus) em uma estrutura hierárquica que divide o espaço recursivamente por eixos alternados (X, Y). Para encontrar o vizinho mais próximo de uma célula, a árvore percorre apenas os ramos relevantes, descartando grandes regiões do espaço sem calcular distâncias desnecessárias.

Exemplo concreto: pontos de ônibus

Os pontos de ônibus geram duas covariáveis: dist_ponto_onibus (distância ao mais próximo) e n_pontos_onibus_250m (contagem em raio de 250m). Ambas capturam a influência do transporte público sobre o ruído — ônibus municipais (frota diesel, sem barreiras acústicas) são fontes sonoras significativas no contexto brasileiro.

Etapa 1: Carregar dados

Os ~11.000 pontos de ônibus são carregados da tabela pontos_onibus (fonte: BHTrans/Prodabel), com coordenadas já em UTM 31983 (metros). A grade de interpolação (~3.3M células) é carregada da tabela grade_interpolacao_10m.

Etapa 2: Construir a árvore

As coordenadas dos pontos de ônibus são inseridas na cKDTree. A construção da árvore leva < 1 segundo para 11K pontos.

Etapa 3: Distância mínima (dist_ponto_onibus)

Para cada célula da grade, a operação tree.query(grade_coords, k=1) encontra o ponto de ônibus mais próximo e retorna a distância euclidiana em metros. A consulta para 3.3 milhões de células leva ~3 segundos.

A célula a 23m de um ponto de ônibus provavelmente sofre impacto acústico direto (frenagem, aceleração, motor diesel em marcha lenta). A célula a 1.847m está em área sem transporte público por ônibus — possivelmente região de preservação ou condomínio fechado.

Etapa 4: Contagem em raio (n_pontos_onibus_250m)

A operação tree.query_ball_point(grade_coords, r=250) conta quantos pontos de ônibus existem dentro de 250m de cada célula. Essa covariável captura a concentração de transporte: um único ponto isolado gera menos ruído que um corredor com 15 pontos em 250m (terminal de integração, eixo viário com faixas exclusivas).

Etapa 5: Transformação e regressão

As distâncias passam por transformação logarítmica (log_e(d + 1)) antes de entrar na regressão. A justificativa é física: a atenuação sonora com a distância segue uma lei logarítmica (decaimento ~6 dB por duplicação de distância em campo livre). A transformação lineariza essa relação, melhorando o ajuste da regressão linear.

As contagens são padronizadas com StandardScaler (média 0, desvio padrão 1) para que covariáveis com escalas diferentes (distância em metros vs. contagem unitária) contribuam de forma equilibrada para o modelo.

Padrão arquitetural: todas as 14 covariáveis

O mesmo padrão cKDTree é aplicado a todas as covariáveis do modelo, mudando apenas o conjunto de pontos de referência:

As duas covariáveis restantes — influencia_via e intensidade_trafego_250m — são compostas: combinam múltiplas distâncias e contagens ponderadas por hierarquia viária e decaimento linear, mas internamente também usam cKDTree.

Performance

O cálculo completo das 14 covariáveis para as ~3.3 milhões de células da grade de BH leva ~45 minutos. A maior parte do tempo é consumida pela intensidade_trafego_250m (soma ponderada com decaimento por distância em raio de 250m) e pela densidade_edif_100m (contagem de ~400K edificações). Sem a cKDTree, o cálculo seria da ordem de dias.

A geração do mapa de ruído municipal envolve processamento pesado sobre ~3.3 milhões de células em grade de 10 m. Esta seção documenta o hardware utilizado e os tempos observados em cada etapa, como referência para reprodutibilidade e planejamento de execuções futuras.

Hardware

Tempos observados por etapa

Todos os tempos referem-se à grade de 10 m (~3.3 M células) e foram medidos nesta máquina.

Tempos totais do pipeline

Considerações de performance

Covariáveis — O gargalo é o cálculo de distâncias e contagens em raio (250 m, 100 m) para cada célula. A estrutura cKDTree (scipy) reduz a complexidade de O(n·m) para O(n·log m), mas o volume — 3.3 M células × centenas de milhares de pontos de referência — ainda demanda ~45 min. O processamento é feito em batches de 20 mil células com commits intermediários, o que permite retomar em caso de interrupção.

Krigagem — A etapa mais custosa. PyKrige opera em single-thread: cada célula resolve um sistema linear com os vizinhos mais próximos. O RK Híbrido adiciona o custo da regressão (Ridge) e do ajuste do variograma sobre os resíduos. A execução de 10 cenários é sequencial e totaliza ~40 h de CPU efetivo.

Rasterização — Etapa leve. Lê os resultados do banco, monta um array NumPy 2D e gera um PNG colormapped com reprojeção EPSG:31983 → 4326 via rasterio.warp. Pode ser reexecutada isoladamente sem recalcular a krigagem.