Zanieczyszczenie powietrza w Polsce - model panelowy
1 Wstęp
Problem zanieczyszczenia powietrza stanowi jedno z najważniejszych wyzwań środowiskowych i zdrowotnych w Polsce. Z mediów od lat płyną ostrzeżenia, że Polska plasuje się w czołówce krajów Unii Europejskiej pod względem stężeń pyłu zawieszonego PM10 - frakcji cząstek stałych o średnicy nieprzekraczającej 10 mikrometrów. Ekspozycja na PM10 wiąże się z istotnym wzrostem ryzyka chorób układu oddechowego i krążenia. Problem ma wyraźny wymiar przestrzenny: stężenia PM10 są znacząco wyższe w miastach przemysłowych południa kraju oraz w regionach o gęstej zabudowie ogrzewanej indywidualnie paliwami stałymi, podczas gdy obszary wiejskie i leśne rejestrują wielokrotnie niższe wartości.
Niniejszy raport analizuje dzienne stężenia PM10 na stacjach pomiarowych GIOŚ w 2024 roku. Główną innowacją metodologiczną jest zastosowanie przestrzennych modeli panelowych, które pozwalają jednocześnie kontrolować stałe cechy stacji, uwzględnić zależności przestrzenne między stacjami oraz modelować dynamikę czasową stężeń. Takie podejście eliminuje błędy wynikające z ignorowania przestrzennej struktury danych.
2 Cel i pytania badawcze
Głównym celem raportu jest identyfikacja i kwantyfikacja determinant dziennych stężeń PM10 na stacjach pomiarowych w Polsce w 2024 roku, z uwzględnieniem przestrzennych zależności między stacjami oraz dynamiki czasowej procesu.
W oparciu o sformułowany cel, postawiłem następujące pytania badawcze:
- W jakim stopniu warunki meteorologiczne wyjaśniają dobową zmienność stężeń PM10, gdy kontroluje się efekty stałe stacji?
- Czy w danych dotyczących jakości powietrza w Polsce występuje istotna autokorelacja przestrzenna, wskazująca na “importowanie” smogu z regionów sąsiednich?
- Która specyfikacja modelu przestrzennego - model błędu przestrzennego (SEM), model opóźnienia przestrzennego (SAR) czy model Cliff-Orda - najlepiej opisuje strukturę zależności przestrzennych w badanych danych?
- Jak silna jest dynamika czasowa stężeń PM10 - na ile stężenie z poprzedniego dnia wyjaśnia stężenie bieżące po kontroli czynników meteorologicznych?
3 Zarys teoretyczny
3.1 Determinanty zanieczyszczenia powietrza
Stężenie pyłów zawieszonych jest wypadkową wielkości lokalnej emisji oraz zdolności atmosfery do jej rozpraszania. Kluczowe zmienne meteorologiczne obejmują:
- Temperaturę powietrza: Wykazuje silnie ujemną korelację z PM10. Spadek temperatury, zwłaszcza w sezonie grzewczym, drastycznie zwiększa zapotrzebowanie na energię cieplną, co napędza tzw. niską emisję z domowych pieców.
- Prędkość wiatru: Działa jako główny mechanizm “czyszczący” (wentylacyjny). Silny wiatr rozprasza zanieczyszczenia i zapobiega ich kumulacji przy powierzchni ziemi.
- Suma opadów: Zjawiska opadowe (deszcz, śnieg) fizycznie “wymywają” cząsteczki PM10 z atmosfery (zjawisko depozycji mokrej).
- Ciśnienie atmosferyczne: Wysokie ciśnienie często wiąże się z bezwietrzną pogodą i osiadaniem mas powietrza, co sprzyja “zamykaniu” zanieczyszczeń w najniższych warstwach atmosfery i powstawaniu zjawiska inwersji temperatury.
Oprócz tego zastosowałem gęstość zaludnienia jako proxy intensywności antropogenicznych źródeł emisji - ruchu drogowego, ogrzewania indywidualnego, działalności przemysłowej.
3.2 Model panelowy
Zastosowanie podejścia panelowego pozwala kontrolować nieobserwowalne cechy stałe w czasie, specyficzne dla każdej stacji - lokalną topografię, dominujące kierunki wiatru, odległość od dróg i zakładów przemysłowych. Efekty stałe stacji absorbują te czynniki, ograniczając ryzyko obciążenia przez pominięte zmienne.
Wybór między modelami dokonywany jest poprzez redukcję od modelu pełnego zgodnie z procedurą Kopczewskiej i in. (2017). Model dynamiczny uwzględniający opóźnienie czasowe PM10 jest uzasadniony fizyką procesu - cząstki zawieszone nie znikają natychmiast po ustaniu emisji, a warunki meteorologiczne zmieniają się powoli.
4 Dane
4.1 Stężenia PM10 - stacje pomiarowe GIOŚ
Dane o stężeniach PM10 pochodzą z bazy pomiarowej Głównego Inspektoratu Ochrony Środowiska (GIOŚ) i obejmują dobowe średnie stężenia pyłu PM10 dla wszystkich aktywnych stacji pomiarowych w Polsce w 2024 roku (\(T\) = 365). Surowe dane zawierały obserwacje dla ponad 160 stacji - jednak ze zbioru wyłączyłem 10 stacji o odsetku braków danych powyżej 10% (\(N\) = 150). Pozostałe braki uzupełniłem interpolacją liniową w czasie (na.approx), ponieważ funkcje pakietu splm wymagają panelu zbilansowanego (Rysunek 1). Finalna baza obejmuje stacje rozmieszczone na terenie całego kraju, przy czym ich zagęszczenie jest silnie niejednorodne - znacznie więcej stacji koncentruje się w miastach i aglomeracjach niż na obszarach wiejskich.
4.2 Gęstość zaludnienia
Dane o zaludnieniu (zmienna pop) pochodzą z Narodowego Spisu Powszechnego 2021 (GUS) w postaci siatki o rozdzielczości 1×1 km. Dla każdej stacji pomiarowej wyznaczyłem bufor o promieniu 5 km, a następnie obliczyłem sumę populacji ze wszystkich komórek siatki zawierających się w tym buforze. Wielkość kropek na Rysunek 2 odzwierciedla wielkość bufora, a ich kolor - populację wewnątrz. Zmienna ta jest stała w czasie (jeden pomiar na stację) i służy jako proxy lokalnej intensywności źródeł antropogenicznych - ruchu drogowego, ogrzewania indywidualnego i działalności przemysłowej. Koniec końców zmienna nie została uwzględniona w modelu, ponieważ wyestymowałem model fixed effects, który ignoruje zmienne stałe w czasie.
4.3 Dane meteorologiczne
Dobowe dane meteorologiczne pobrałem z serwisu Open-Meteo, który udostępnia dane bardzo dużej rozdzielczości (1 - 11 km). Dzięki temu uzyskałem dane osobno dla współrzędnych każdej stacji pomiarowej. Pobrałem pięć zmiennych dobowych za okres 1 stycznia - 31 grudnia 2024:
temperature_2m_mean: średnią temperaturę powietrza na wysokości 2 m (°C),precipitation_sum: sumę opadów (mm),wind_speed_10m_max: maksymalną prędkość wiatru na wysokości 10 m (km/h),surface_pressure_mean: średnie ciśnienie atmosferyczne przy powierzchni (hPa),relative_humidity_2m_mean: średnią wilgotność względną na wysokości 2 m (%)
Przypisanie danych meteorologicznych do konkretnych współrzędnych stacji (zamiast korzystania z danych dla najbliższej stacji synoptycznej) eliminuje błąd przestrzenny wynikający z odległości między stacją PM10 a stacją meteorologiczną. W danych meteorologicznych nie stwierdziłem braków.
Przed przystąpieniem do modelowania sprawdziłem, czy w danych występuje trend czasowy. Charakterystyka modeli panelowych (oszacowanie jednego współczynnika w czasie) wymaga, aby trend czasowy był stabilny i jednokierunkowy. Rysunek 3 przedstawia przebieg w czasie zmiennych wykorzystanych w modelu z uwzględnieniem rozstępu między kwantylami 0.05 i 0.95. Zmienna modelowana pm10 wykazuje większą wariancję w sezonie “grzewczym”, natomiast nie widać znaczącego tredu. Nieliniowy trend występuje w przypadku zmiennych temperature_2m_mean i relative_humidity_2m_mean, co można powiązać z cyklem pór roku.
5 Modelowanie
5.1 Macierz wag przestrzennych
Jak widać na mapie (Rysunek 2), stacje pomiarowe nie są rozłożone równomiernie w kraju. Duża część z nich jest skoncentrowana wokół dużych ośrodków miejskich, natomiast obszary mniej zurbanizowane - jak Podlasie i Mazury - mają tylko kilka stacji. Ta obserwacja wpływa na wybór odpowiedniej macierzy wag przestrzennych. Na przykład macierz \(k\) najbliższych sąsiadów zniekształciłaby rzeczywiste zależności przestrzenne, ponieważ najbliżsi sąsiedzi odizolowanej stacji mieliby taką samą wagę jak stacje położone blisko siebie w dużym mieście. Z tego względu zastosowałem macierz odwrotnej odległości.
Balanced Panel: n = 150, T = 365, N = 54750stations_ordered <- stations_pm10_df$gios_code
stations_sf <- st_as_sf(stations_pm10_df,
coords = c("lon", "lat"),
crs = 4326) %>%
st_transform(crs = 2180)
dist_mx <- st_distance(stations_sf)
dist_mx <- as.numeric(dist_mx) / 1000
dist_mx <- matrix(dist_mx,
nrow = nrow(stations_pm10_df),
ncol = nrow(stations_pm10_df))
w_mx <- 1 / (dist_mx^1 + 1) # dodanie 1 w mianowniku dla uniknięcia błędów
diag(w_mx) <- 0
w_mx <- mat2listw(w_mx, style = "W")5.2 Diagnostyka
Formę funkcyjną modelu opisuje równanie eq_1:
eq_1 <- pm10 ~ temperature_2m_mean +
pm10_lag1 +
precipitation_sum +
surface_pressure_mean +
relative_humidity_2m_mean +
wind_speed_10m_maxNiestety pomimo poprawnej specyfikacji macierzy wag przestrzennych i panelowej ramki danych pdata.frame, testy diagnostyczne bsjktest, bsktest i sphtest z pakietu splm nie zadziałały. Błąd polega na niezgodności wymiarów mnożonych macierzy wewnątrz funkcji. Debuggowanie nie przyniosło skutku, więc przeszedłem do kolejnego etapu bez szczegółowych testów.
tic("LMH"); bsktest(eq_1, data = data_pdf, listw = w_mx, test = "LMH"); toc()Error in `WWp %*% q`:
! non-conformable argumentsJedyny działający test LM1 pozwala stwierdzić, że efekty losowe są obecne i nie należy stosować regresji łącznej (pooled regression).
tic("LM1"); bsktest(eq_1, data = data_pdf, listw = w_mx, test = "LM1"); toc()
Baltagi, Song and Koh SLM1 marginal test
data: pm10 ~ temperature_2m_mean + pm10_lag1 + precipitation_sum + surface_pressure_mean + relative_humidity_2m_mean + wind_speed_10m_max
LM1 = 851.22, p-value < 2.2e-16
alternative hypothesis: Random effectsLM1: 1.67 sec elapsedBez znajomości statystyk testowych uznałem, że bezpieczniejszym wyborem będzie zastosowanie estymatora fixed effects. Model FE jest zazwyczaj stosowany gdy wymiar \(N\) panelu jest niewielki, a w tym przypadku tak jest (Kopczewska, Kudła, i Walczyk 2017).
5.3 Wybór modelu
Zgodnie z przyjętą procedurą, rozpocząłem od wyestymowania pełnego modelu Cliffa-Orda fixed effects z opóźnieniem przestrzennym zmiennej modelowanej i błędu. Dodatkowo do modelu włączyłem opóźnioną czasowo zmienną pm10_lag1, która pozwala uchwycić efekt utrzymywania się smogu przez kolejne dni. Taki zabieg w modelach FE prowadzi do obciążenia Nickella (Nickell 1981), wynikającego z faktu, że estymator FE powstaje przez odjęcie średniej. Jednak przy dużym wymiarze czasowym (T = 365), obciążenie jest zaniedbywalne.
W następnym kroku wyestymowałem modele z nałożonymi ograniczeniami na komponenty przestrzenne: SAR i SEM. Niestety nie mogłem porównać ich za pomocą kryteriów AIC, ponieważ model Cliff-Orda w pakiecie splm jest estymowany metodą quasi-największej wiarygodności i obiekt modelu nie zawiera wartości funkcji wiarygodności logLik.
Z analizy wydruku modelu Cliff-Orda możemy stwierdzić, że oszacowania współczynników przestrzennych (\(\rho\) i \(\lambda\)) są silnie istotne. Co więcej, oba współczynniki są dodatnie, co jest zgodne z intuicją i nie budzi podejrzeń dotyczących \(overfittingu\) modelu. Zdecydowałem się zatem przyjąć model Cliff-Orda za najbardziej odpowiedni.
model_sarar <- spml(eq_1,
data = data_pdf,
listw = w_mx,
model = "within",
effect = "individual",
lag = TRUE,
spatial.error = "b")
summary(model_sarar)Spatial panel fixed effects sarar model
Call:
spml(formula = eq_1, data = data_pdf, listw = w_mx, model = "within",
effect = "individual", lag = TRUE, spatial.error = "b")
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-97.6067 -3.1293 0.7184 3.9223 159.8582
Spatial error parameter:
Estimate Std. Error t-value Pr(>|t|)
rho 0.9238875 0.0052426 176.23 < 2.2e-16 ***
Spatial autoregressive coefficient:
Estimate Std. Error t-value Pr(>|t|)
lambda 0.802855 0.011214 71.597 < 2.2e-16 ***
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
temperature_2m_mean -0.2827351 0.0290271 -9.7404 < 2.2e-16 ***
pm10_lag1 0.4556932 0.0036016 126.5252 < 2.2e-16 ***
precipitation_sum -0.0469780 0.0115051 -4.0832 4.441e-05 ***
surface_pressure_mean -0.1786357 0.0206206 -8.6630 < 2.2e-16 ***
relative_humidity_2m_mean -0.0857089 0.0073970 -11.5869 < 2.2e-16 ***
wind_speed_10m_max -0.2100764 0.0098160 -21.4014 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.4 Efekty indywidualne
Jak widać na poniższych wykresach, rozkład efektów indywidualnych przypomina krzywą gaussowską, jednak ma grubszy lewy ogon. Mapa efektów indywidualnych pokazuje ciekawą zależność. Stacje położone bardziej na południe (a więc na wyżynach i w górach), mają niższe wartości PM10 niż przewiduje model. Efekt indywidualny reprezentuje wpływ czynników niezmiennych w czasie, jak również czynników pominiętych w modelu. W dużej mierze są to czynniki antropogeniczne. Ich wyjaśnienie stanowi przyczynek do pogłębienia tematu.
eff <- effects(model_sarar)
ind_eff <- eff$SETable[, 1]
names(ind_eff) <- NULL
stations_sf$ind_eff <- ind_eff
plot(density(ind_eff), main = "Rozkład efektów indywidualnych")
ggplot() +
geom_sf(data = woj_sf) +
geom_sf(data = stations_sf, aes(color = ind_eff)) +
scale_color_gradientn(colors = c("#1E88E5", "#FFEB3B", "#D32F2F")) +
theme_minimal() +
labs(color = "Efekt indywidualny",
title = "Mapa efektów indywidualnych")
5.5 Efekty krańcowe
Zgodnie z teorią ekonometrii przestrzennej, w sytuacjach, gdy w modelu uwzględniono przestrzenne opóźnienie zmiennej zależnej (tzw. efekt równoczesności), wyników nie należy interpretować wprost za pomocą standardowych współczynników \(\beta\), lecz z wykorzystaniem miar wpływu: efektów bezpośrednich (direct) oraz pośrednich (indirect). Wyliczenie tych miar jest konieczne, aby uchwycić zjawisko lokalnej wariancji w czasie i przestrzeni.
attr(model_sarar, "have_factor_preds") <- FALSE
impacts(model_sarar, listw = w_mx, time = 10, R = 1000)Impact measures (lag, trace):
Direct Indirect Total
temperature_2m_mean dy/dx -0.28299753 -1.1491781 -1.4321757
pm10_lag1 dy/dx 0.45611616 1.8521671 2.3082833
precipitation_sum dy/dx -0.04702163 -0.1909424 -0.2379640
surface_pressure_mean dy/dx -0.17880154 -0.7260658 -0.9048673
relative_humidity_2m_mean dy/dx -0.08578845 -0.3483642 -0.4341527
wind_speed_10m_max dy/dx -0.21027140 -0.8538566 -1.0641280Dla każdej zmiennej objaśniającej w modelu efekty pośrednie są średnio 4-krotnie wyższe od efektów bezpośrednich, co pokazuje skalę oddziaływań przestrzennych. Stężenie smogu z dnia poprzedniego jest najsilniejszym stymulantem dzisiejszego zanieczyszczenia (Total = 2.308). Wskazuje to na zjawisko kumulacji zanieczyszczeń w atmosferze. Temperatura powietrza wywiera najsilniejszy ujemny wpływ na PM10 spośród wszystkich zmiennych pogodowych (Total = -1.432). Zgodnie z teorią, spadek temperatury napędza tzw. niską emisję (zwiększone spalanie paliw stałych w sektorze komunalno-bytowym), co drastycznie pogarsza jakość powietrza lokalnie i regionalnie. Silny wiatr działa jako czynnik wentylujący, skutecznie obniżając stężenie pyłów (Total = -1.064).
Spośród oszacowań efektów krańcowych, tylko współczynnik dla ciśnienia atmosferycznego odbiega od przewidywań teoretycznych. Zimowe wyże baryczne sprzyjają zamykaniu smogu przy ziemi (spodziewany znak dodatni). Otrzymane oszacowanie wymaga dyskusji i sygnalizuje potrzebę dalszego udoskonalenia modelu.
6 Podsumowanie
Wyniki wskazują, że warunki meteorologiczne w istotnym stopniu kształtują poziom zanieczyszczeń - wyższa temperatura, silniejszy wiatr oraz opady atmosferyczne sprzyjają spadkowi stężenia pyłu. Najsilniejszym czynnikiem wyjaśniającym bieżące stężenie PM10 jest jednak jego poziom z dnia poprzedniego, co potwierdza dużą trwałość zanieczyszczeń w atmosferze.
Analiza wykazała również wyraźną autokorelację przestrzenną, co oznacza, że zanieczyszczenia są transportowane między regionami i poziom PM10 w jednej lokalizacji zależy od sytuacji w lokalizacjach sąsiednich. Spośród rozważanych specyfikacji najlepiej strukturę danych opisuje model Cliffa–Orda, który uwzględnia zarówno przestrzenne oddziaływania między stacjami, jak i przestrzenną zależność składnika losowego. W pracy zastosowałem macierz wag przestrzennych opartą na odwrotnej odległości między stacjami, co pozwala lepiej odwzorować rzeczywisty proces transportu zanieczyszczeń w sytuacji nierównomiernego rozmieszczenia stacji pomiarowych.