---
title: "Drooni Uudellamaalla: Paljonko hälytys maksoi suhteessa todelliseen riskiin?"
subtitle: "Bayesilainen ROI-analyysi 15.5.2026 vaaratiedotteesta"
date: "2026-05-15"
author: "Kristian Vepsäläinen"
categories: [riskianalyysi, bayesilainen tilastotiede, avoin data, turvallisuus]
code-fold: true
execute:
warning: false
message: false
format:
html:
theme: minimal
fig-width: 9
fig-height: 5.5
toc: true
---
```{r setup}
library(tidyverse)
library(fst)
library(here)
library(spdep)
library(brms)
library(scales)
library(patchwork)
# Väripaletti
COL_RED <- "#e63946"
COL_GREEN <- "#2a9d8f"
COL_ORANGE <- "#f4a261"
COL_NAVY <- "#1d3557"
COL_BLUE <- "#457b9d"
```
## Tänä aamuna Suomi heräsi droonihälytykseen
Perjantaina 15.5.2026 kello 4 yöllä Uudellamaalla annettiin vaaratiedote: Puolustusvoimien uhka-arvion perusteella alueelle saattoi suuntautua vaarallinen drooni. Helsinki-Vantaan lentoliikenne pysäytettiin, Hornetteja lensi Helsingin yllä, ja tuhansia ihmisiä kehotettiin hakeutumaan sisätiloihin.
Kysymys ei ole siitä, toimivatko viranomaiset oikein — ehkä toimivat. Kysymys on: **mitä tämä oikeasti maksoi, ja mihin riskiin suhteutettuna?**
Tämä on se analyysi, jota kukaan muu ei tee, koska se vaatii sekä tilastollista osaamista että kykyä hakea avoimesta datasta oikeat luvut. Pistearvo ("drooni havaittiin") on arvoton ilman jakaumaa. **Maailma on jakauma.**
---
## 1. Missä tämä drooni olisi voinut räjähtää? — Paavo-data Uudeltamaalta
Data on ladattu Tilastokeskuksen Paavo-tietokannasta (2024) ja tallennettu paikallisesti `.fst`-muotoon.
```{r paavo-data}
# Lue paikallisesti tallennettu Paavo-data
paavo_raw <- read_fst(here("data", "paavo", "paavo.fst"))
df_uusimaa <- paavo_raw |>
filter(Vuosi == 2024) |>
# Uusimaa: postinumerot alkavat 0 (00xxx–09xxx) tai 10 (Hyvinkää, Järvenpää yms.)
filter(str_starts(Postinumeroalue, "0") | str_starts(Postinumeroalue, "10")) |>
as_tibble() |>
rename(
postinumero = Postinumeroalue,
vaesto = `Asukkaat yhteensä (HE)`,
asunnot = `Asunnot (RA)`,
tyopaikat = `Työpaikat yhteensä (TP)`,
pinta_ala = `Postinumeroalueen pinta-ala`
) |>
mutate(across(c(vaesto, asunnot, tyopaikat, pinta_ala), as.numeric)) |>
filter(!is.na(vaesto), vaesto > 0) |>
mutate(
# Paavo ilmoittaa pinta-alan neliökilometreinä (km²)
# Diagnostiikka alla varmistaa suuruusluokan
pinta_ala_km2 =pinta_ala/1e6
)
# Diagnostiikka — varmistaa että pinta_ala on km² eikä m²
# Uudenmaan oikea pinta-ala ~9 400 km²; jos summa ~9 400e6 → yksikkö on m²
stopifnot(
"pinta_ala näyttää olevan m² eikä km² — tarkista yksikkö" =
sum(df_uusimaa$pinta_ala_km2, na.rm = TRUE) < 1e5
)
cat("pinta_ala_km2 summa:", round(sum(df_uusimaa$pinta_ala_km2, na.rm = TRUE)), "km²",
"(Uudenmaan oikea ~9 400 km²)\n")
cat("Postinumeroalueita:", nrow(df_uusimaa), "\n")
cat("Väestö yhteensä: ", format(sum(df_uusimaa$vaesto, na.rm = TRUE), big.mark = " "), "\n")
cat("Asuntoja yhteensä: ", format(sum(df_uusimaa$asunnot, na.rm = TRUE), big.mark = " "), "\n")
```
```{r uusimaa-vis}
# Väestöjakauma postinumeroalueittain
p1 <- df_uusimaa |>
arrange(desc(vaesto)) |>
slice_head(n = 30) |>
mutate(postinumero = fct_reorder(postinumero, vaesto)) |>
ggplot(aes(x = vaesto, y = postinumero)) +
geom_col(fill = COL_NAVY, alpha = 0.85) +
scale_x_continuous(labels = label_number(big.mark = " ")) +
labs(
title = "30 väkirikkainta postinumeroaluetta Uudellamaalla",
x = "Asukkaita", y = NULL,
caption = "Lähde: Tilastokeskus Paavo 2023"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
p1
```
---
## 2. Todennäköisyys: mihin satunnainen drooni osuu?
Tarkastellaan tilannetta kahdella tasolla: **geometrinen malli** (pinta-ala) ja **väestöpainotettu malli** (ihmisiä per neliökilometri).
```{r drone-probability}
# Uudenmaan kokonaispinta-ala Paavo-datasta
pinta_ala_yht <- sum(df_uusimaa$pinta_ala_km2, na.rm = TRUE) # km² (yksikkökorjattu)
vaesto_yht <- sum(df_uusimaa$vaesto, na.rm = TRUE)
# Droonin "osumisala": räjähdyssäde perustuu dokumentoituihin tietoihin
#
# Shahed-136 / Geran-2 -tyyppi (venäläis-iranilainen, laajasti käytetty Ukrainassa):
# - Taistelukärki: 40–50 kg räjähdysainetta (Wikipedia: HESA Shahed 136;
# Army Recognition; ISIS-Online Alabuga-raportti 2024)
# - Räjähdyspaineaalto: merkittäviä rakennevaurioita 15–30 m säteellä
# - Sirpaleiden tappava säde: 30–50 m, vammauttava säde yli 100 m kaupunkialueella
# (Lähde: shahed136.com tekninenanalyysi; Wikipedia fragmentaatiokuvaus)
# - Suora osuma rakennukseen: kymmeniä uhreja mahdollinen (rakenteen romahdus,
# tulipalo, fragmentaatio) — ISIS-Online: "direct hits on occupied buildings
# can cause dozens of casualties"
#
# Käytetään konservatiivista tappavaa sädettä 50 m (sirpaleuhka),
# joka on lähteen alaraja-arvio. Yläarvio (100 m) käytetään herkkyysanalyysissä.
r_m <- 50 # räjähdyssäde metreinä — konservatiivinen alaraja
r_km <- r_m / 1000
osuma_ala_km2 <- pi * r_km^2
# P(osuu mihinkään rakennukseen) — geometrinen malli
# Rakennettu pinta-ala Uudellamaalla: ~15–20% kokonaispinta-alasta (kaupunkialue)
# Käytetään konservatiivista arviota 20%
rakennettu_osuus <- 0.20
rakennettu_ala <- pinta_ala_yht * rakennettu_osuus
# Geometrinen todennäköisyys: drooni osuu johonkin 50m säteellä
p_geometrinen <- osuma_ala_km2 / pinta_ala_yht
# Väestötiheys Uudellamaalla
tiheys_km2 <- vaesto_yht / pinta_ala_yht # henkilöä/km²
# ── Kaksi eri kysymystä, kaksi eri suuretta ────────────────────────────────
#
# KYSYMYS A: E[kuolemia | drooni räjähtää jossakin Uudellamaalla]
# = tiheys × osuma_ala × rakennettu_osuus × p_kuolema_per_hlö
# Tulkinta: odotusarvo kuolemien lukumäärästä.
# p_geometrinen EI kuulu tähän — drooni JO räjähti jossakin,
# ja kysymme kuinka monta kuollee keskimäärin.
#
# KYSYMYS B: P(≥1 kuolema | drooni räjähtää jossakin Uudellamaalla)
# Poisson-approksimaatio: P(X≥1) = 1 - exp(-lambda), lambda = E[kuolemia]
# Koska lambda << 1, tämä ≈ lambda numeerisesti — mutta käsitteellisesti eri.
#
# Virhe aiemmin: p_geometrinen × ihmiset_sateessa kertolasku sisälsi
# osuma_ala kahdesti (kerran p_geo:ssa, kerran ihmiset_sateessa:ssa).
# ───────────────────────────────────────────────────────────────────────────
p_kuolema_per_henkilo <- 0.30 # P(kuolema | ihminen säteessä, yöaika, sisätilat)
# A: odotusarvo kuolemia
ev_kuolemia <- tiheys_km2 * osuma_ala_km2 * rakennettu_osuus * p_kuolema_per_henkilo
# B: todennäköisyys vähintään yhdelle kuolemalle
p_vah_1_kuolema <- 1 - exp(-ev_kuolemia)
tibble(
Skenaario = c(
"P(drooni osuu 50m säteellä johonkin kohtaan Uudellamaalla)",
"P(drooni osuu asuinalueelle) [geo. × rakennettu 20%]",
"E[ihmisiä räjähdyssäteessä] [tiheys × ala × rakennettu]",
"E[kuolemia | drooni räjähtää jossakin] — Kysymys A",
"P(≥1 kuolema | drooni räjähtää jossakin) — Kysymys B"
),
Arvo = c(
sprintf("%.2e", p_geometrinen),
sprintf("%.2e", p_geometrinen * rakennettu_osuus),
sprintf("%.4f hlö", tiheys_km2 * osuma_ala_km2 * rakennettu_osuus),
sprintf("%.4f kuollutta", ev_kuolemia),
sprintf("%.4f (%.2f%%)", p_vah_1_kuolema, p_vah_1_kuolema * 100)
)
) |>
knitr::kable(caption = "Todennäköisyysarviot — drooni räjähtää Uudellamaalla")
```
### Bayesilainen analyysi: epävarmuus huomioituna
Pelkät pisteet ovat harhaanjohtavia. Käytetään Monte Carlo -simulaatiota, joka ottaa huomioon epävarmuuden kaikissa parametreissa.
```{r bayesian-monte-carlo}
set.seed(20260515)
N <- 100000
# Epävarmat parametrit jakaumina
sim <- tibble(
# Räjähdyssäde: dokumentoitu vaihteluväli Shahed-tyyppisille drooneille
# 30–50 m tappava sirpalevyöhyke, 50–100 m vammauttava vyöhyke
# (shahed136.com; ISIS-Online; Army Technology)
r_m_sim = runif(N, 30, 100),
# Rakennettu osuus: 15–35%
rak_osuus_sim = rbeta(N, 4, 12) + 0.10,
# Tappavuus yöaikaan ulkona: 0.05–0.40 (vähemmän ihmisiä ulkona klo 4)
tappavuus_sim = rbeta(N, 2, 8),
# Väestötiheys: epävarmuus ±20%
tiheys_sim = rnorm(N, tiheys_km2, tiheys_km2 * 0.15)
) |>
mutate(
ala_sim = pi * (r_m_sim / 1000)^2,
# p_geo_sim: geometrinen osuma — EI käytetä E[kuolemia]-laskussa
p_geo_sim = ala_sim / pinta_ala_yht,
# ihmiset: E[ihmisiä säteessä | räjähtää jossakin]
ihmiset = pmax(0, tiheys_sim * ala_sim * rak_osuus_sim),
# ev_kuolemia: E[kuolemia | räjähtää jossakin] — p_geo_sim EI mukana
# (se lisättäisiin vain jos laskettaisiin E yli kaikkien mahdollisten räjähdyspaikkojen)
ev_kuolemia = ihmiset * tappavuus_sim,
# P(≥1 kuolema) Poisson-approksimaatiolla
p_kuolema = 1 - exp(-ev_kuolemia)
)
# Posteriorijakauma
p1_mc <- sim |>
ggplot(aes(x = p_kuolema)) +
geom_histogram(bins = 80, fill = COL_RED, alpha = 0.75) +
geom_vline(
xintercept = quantile(sim$p_kuolema, c(0.05, 0.50, 0.95)),
linetype = c("dashed","solid","dashed"),
color = COL_NAVY, linewidth = 0.8
) +
scale_x_log10(labels = label_scientific()) +
labs(
title = "Posteriorijakauma: E[kuolemia | drooni räjähtää satunnaisesti Uudellamaalla]",
subtitle = sprintf(
"Mediaani: %.2e | 90%% credible interval: [%.2e, %.2e]",
median(sim$p_kuolema),
quantile(sim$p_kuolema, 0.05),
quantile(sim$p_kuolema, 0.95)
),
x = "Todennäköisyys (log-asteikko)",
y = "Simulaatioita"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
p1_mc
```
---
## 3. Vertailu: työmatkaliikenne
Kuinka tämä riski suhteutuu johonkin, minkä kanssa elämme päivittäin?
```{r tyomatkariski}
# Lähde: Liikenneturva / Tilastokeskus
# Tieliikennekuolemia Suomessa ~200/vuosi
# Työmatkakuolemia n. 10–15% kaikista = ~20–30/vuosi
# Työssäkäyviä Uudellamaalla ~800 000 (Paavo-data)
# Työmatkoja per vuosi per henkilö: ~240 (5 pv/vko × 48 vko)
# Työmatkakuolemia per suoritettu matka
tyomatka_kuolemat_v <- 25 # Suomi/vuosi, konservatiivinen
tyossaakayvia_suomi <- 2.4e6 # työllistettyä
matkat_v <- tyossaakayvia_suomi * 2 * 240 # meno + paluu
p_kuolema_tyomatka_per_matka <- tyomatka_kuolemat_v / matkat_v
# Uudenmaan osuus ~1/3 Suomen työssäkäyvistä
# Verrataan: 1 drooni vs 1 työmatkavuosi (2 matkaa/päivä, 240 päivää)
matkat_vuodessa <- 480
riski_tyomatka_vuosi <- 1 - (1 - p_kuolema_tyomatka_per_matka)^matkat_vuodessa
# Vertailutaulukko käyttää kahta eri suuretta:
# ev_kuolemia = E[kuolemia | räjähtää jossakin] — odotusarvo kuolemista
# p_vah_1_kuolema = P(≥1 kuolema | räjähtää jossakin) — todennäköisyys
# Työmatkariski on yksilön vuosiriski — lähimpänä p_vah_1_kuolema:a.
tibble(
Tilanne = c(
"Drooni räjähtää Uudellamaalla: E[kuolemia]",
"Drooni räjähtää Uudellamaalla: P(≥1 kuolema)",
"Drooni räjähtää: P(≥1 kuolema), MC mediaani",
"Kuolet yhden työmatkavuoden aikana (Suomi)",
"Kuolet tupakoinnin takia (tupakoitsija, 20v)",
"Kuolet salamaniskuun Suomessa (vuosi)"
),
Arvo = c(
ev_kuolemia,
p_vah_1_kuolema,
median(sim$p_kuolema),
riski_tyomatka_vuosi,
0.015,
1/1e6
),
Yksikkö = c(
"E[kuolemia per räjähdys]",
"P(≥1 kuolema per räjähdys)",
"P(≥1 kuolema per räjähdys, MC)",
"P(kuolema / hlö / vuosi)",
"P(kuolema / hlö / vuosi)",
"P(kuolema / hlö / vuosi)"
)
) |>
arrange(desc(Arvo)) |>
mutate(Arvo = sprintf("%.2e", Arvo)) |>
knitr::kable(caption = "Riskivertailu: drooni vs. arkipäivän riskit")
```
```{r riskivertailu-vis}
risks <- tibble(
label = c(
"Työmatkavuosi\n(liikennekuolema)",
"Drooni räjähtää\n(95. persentiili)",
"Drooni räjähtää\n(mediaani)",
"Salamanisku\n(1 vuosi)"
),
p = c(
riski_tyomatka_vuosi,
as.numeric(quantile(sim$p_kuolema, 0.95)),
median(sim$p_kuolema),
1e-6
),
kategoria = c("Taustariskit", "Drooni", "Drooni", "Taustariskit")
)
risks |>
mutate(label = fct_reorder(label, p)) |>
ggplot(aes(x = p, y = label, fill = kategoria)) +
geom_col(alpha = 0.85) +
geom_label(aes(label = sprintf("%.1e", p)), hjust = -0.1, size = 3.5, fill = "white") +
scale_x_log10(labels = label_scientific(), limits = c(1e-9, 1)) +
scale_fill_manual(values = c("Drooni" = COL_RED, "Taustariskit" = COL_BLUE)) +
labs(
title = "Riski kuolla: drooni vs. arkipäivän vaarat",
subtitle = "Drooni-riski on häiriöskenaariossa pienempi kuin työmatkavuosi",
x = "P(kuolema) — log-asteikko",
y = NULL,
fill = NULL,
caption = "Lähteet: Liikenneturva, Tilastokeskus, oma simulaatio"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top"
)
```
---
## 4. ROI-analyysi: Paljonko hälytys maksoi?
Nyt lasketaan toimenpiteen kustannukset suhteessa estettyyn vahinkoon.
```{r roi-laskenta}
# ---- KUSTANNUSPUOLI ----
# Helsinki-Vantaa suljettu n. klo 04:00–08:00 = ~4 tuntia
# 2024: 16,3 milj. matkustajaa → ~44 600/pv → ~1 860/h
# Lennon viivästyminen/peruutus: matkustaja menettää aikaa + lentoyhtiöt saavat
# korvausvastuita (EU 261/2004: 250–600 € / matkustaja peruttaessa)
hki_matkustajia_pv <- 16.3e6 / 365
hki_matkustajia_h <- hki_matkustajia_pv / 24
sulkemisaika_h <- 4
hairiintyneet_matk <- hki_matkustajia_h * sulkemisaika_h
# Välittömät kustannukset
lentoyhtion_tappio_per_matk <- 180 # € (polttoaine, henkilöstö, gate-maksut, kompensaatiot)
matkustajan_aika_per_matk <- 45 # € (2h viive × 22,5 €/h, Suomen BOT-arvio)
kust_lentokentta <- hairiintyneet_matk * (lentoyhtion_tappio_per_matk + matkustajan_aika_per_matk)
# HUS: leikkauksia siirretty / henkilökuntaa ei saapunut ajoissa
# HUSin budjetti 2024 ~2,4 mrd → ~6,6 milj./vrk → ~275 000 €/h
# Viivästyneet leikkaukset klo 5–9: n. 4h kriittistä aamutoimintaa
# Konservatiivinen arvio: 10% aamun kapasiteetista häiriintyi = 2h efektiivistä hukan
hus_vrk_budjetti <- 2.4e9 / 365
hus_h_budjetti <- hus_vrk_budjetti / 24
kust_hus <- hus_h_budjetti * 2 * 0.10 # 2h, 10% kapasiteetti
# Puolustusvoimat: Hornet-lento
# F/A-18C Hornet lentotunti ~10 000 € (NATO-arviot)
hornet_lentotunteja <- 4 * 2 # 2 konetta, 4h
kust_hornetit <- hornet_lentotunteja * 10000
# Pelastuslaitosten ylitöitä, viranomaisviestintä yms.
kust_muut <- 200000
kust_yhteensa <- kust_lentokentta + kust_hus + kust_hornetit + kust_muut
# ---- HYÖTYPUOLI ----
# Estetty vahinko, JOS drooni olisi räjähtänyt keskeisessä kohteessa
# Worst case: ihmishengenmenetykset
# Value of Statistical Life (VSL) Suomessa: ~3–4 milj. € (Väylävirasto 2021)
vsl_eur <- 3.5e6
# Odotettu kuolleisuus ilman toimenpidettä:
# p(drooni on todellinen uhka) × p(räjähdys) × p(kuolema)
# — viranomaiset eivät tiedä, onko drooni todellinen vai vääräpositiivi
# — historiallisesti vastaavista hälytyksistä: n. 5–15% johtaa todelliseen tilanteeseen
# Käytetään optimistista uhka-arviota: p(todellinen) = 0.15
p_todellinen_uhka <- 0.15
odotetut_kuolemat_ilman <- p_todellinen_uhka * 0.3 * 2 # 2 mahdollista uhria (keskim.)
estetty_vahinko <- odotetut_kuolemat_ilman * vsl_eur
# ROI = (Estetty vahinko - Kustannukset) / Kustannukset
roi <- (estetty_vahinko - kust_yhteensa) / kust_yhteensa
tibble(
Erä = c(
"Lentokenttä (matkustajat + lentoyhtiöt)",
"HUS-kapasiteettihäiriö",
"Hornet-lentotunnit",
"Muut viranomaiskulut",
"KUSTANNUKSET YHTEENSÄ",
"",
"Odotettu estetty vahinko (p=0.15 todellinen uhka)",
"NETTOHYÖTY",
"ROI"
),
`Euroa (€)` = c(
kust_lentokentta,
kust_hus,
kust_hornetit,
kust_muut,
kust_yhteensa,
NA,
estetty_vahinko,
estetty_vahinko - kust_yhteensa,
NA
),
Huomio = c(
sprintf("%d matkustajaa × %d €/matk.", round(hairiintyneet_matk), lentoyhtion_tappio_per_matk + matkustajan_aika_per_matk),
"10% 2h kapasiteetti",
sprintf("%d lentotuntia × 10k€", hornet_lentotunteja),
"Arvio",
"",
"",
sprintf("p=%.2f × 2 henkeä × %.1f M€/VSL", p_todellinen_uhka, vsl_eur/1e6),
"",
sprintf("%.2f (negatiivinen = häviöllinen)", roi)
)
) |>
filter(!is.na(`Euroa (€)`) | Erä %in% c("KUSTANNUKSET YHTEENSÄ", "NETTOHYÖTY")) |>
mutate(`Euroa (€)` = case_when(
Erä == "ROI" ~ NA_real_,
TRUE ~ `Euroa (€)`
)) |>
mutate(`Euroa (€)` = format(round(`Euroa (€)`), big.mark = " ", nsmall = 0)) |>
knitr::kable(caption = "ROI-laskelma: Droonivastatoimet 15.5.2026")
```
### ROI:n herkkyysanalyysi
ROI riippuu kriittisesti kahdesta parametrista: **uhkan aitouden todennäköisyydestä** ja **odotettujen uhrien määrästä**.
```{r roi-herkkyys}
p_uhka_grid <- seq(0.01, 0.50, length.out = 50)
uhrit_grid <- c(1, 3, 10, 30)
roi_df <- expand_grid(p_uhka = p_uhka_grid, uhrit = uhrit_grid) |>
mutate(
estetty = p_uhka * 0.3 * uhrit * vsl_eur,
roi_val = (estetty - kust_yhteensa) / kust_yhteensa
)
roi_df |>
mutate(uhrit = factor(uhrit, labels = paste0(uhrit_grid, " uhria"))) |>
ggplot(aes(x = p_uhka, y = roi_val, color = uhrit)) +
geom_line(linewidth = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
annotate("rect",
xmin = 0.01, xmax = 0.50,
ymin = -Inf, ymax = 0,
fill = COL_RED, alpha = 0.05
) +
annotate("text", x = 0.48, y = -0.7, label = "Häviöllinen", color = COL_RED, hjust=1) +
annotate("text", x = 0.48, y = 0.3, label = "Kannattava", color = COL_GREEN, hjust=1) +
scale_x_continuous(labels = percent_format()) +
scale_color_manual(values = c(COL_BLUE, COL_GREEN, COL_ORANGE, COL_RED)) +
labs(
title = "ROI-herkkyysanalyysi: Milloin hälytys on taloudellisesti perusteltu?",
subtitle = sprintf("Kiinteät kustannukset: %.1f M€ | VSL = %.1f M€", kust_yhteensa/1e6, vsl_eur/1e6),
x = "P(uhka on todellinen)",
y = "ROI",
color = "Odotetut uhrit",
caption = "Lähteet: Finavia 2024, HUS, Väylävirasto VSL, oma laskelma"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "right"
)
```
---
## 4b. Fat tail: entä jos drooni pudottaa täyteen tankatun matkustajakoneen kaupungin päälle?
Perusskenaario käsitteli droonia, joka räjähtää satunnaisessa paikassa. Mutta on olemassa toinen, paljon vakavampi skenaario: **drooni osuu lähestyvään matkustajakoneeseen matalalentovaiheessa polttoaine täynnä**. Tämä on juuri se tilanne, jonka ilmatilansulku estää suoraan.
Ennen kuin lasketaan, on syytä tarkistaa geometria: menevätkö lähestymisreitit Helsingin kaupungin päälle?
Helsinki-Vantaalla on kolme kiitotietä. Ensisijaisesti laskeudutaan kiitotielle 2 (15) luoteen eli Nurmijärven suunnasta tai kiitotielle 1 (22L) koillisesta Keravan suunnasta. Kiitotienumero kertoo kompassisuunnan: kiitotie 15 = 150°, kiitotie 22L = 220°. Lähestyminen tapahtuu vastakkaisesta suunnasta — kiitotielle 15 tullaan siis koillisesta (suunta ~330°, Nurmijärven suunnasta), kiitotielle 22L tullaan lounaasta (suunta ~40°, Keravan suunnalta).
Kumpikaan reitti ei lennä suoraan Helsingin kantakaupungin yli. Vantaan kaupunkialue sen sijaan osuu molempiin lähestymiskäytäviin. Tämä on tärkeä havainto: **worst-case-skenaario ei ole Helsingin keskusta vaan Vantaan tiiviit asuinalueet** lähestymiskäytävien alla.
```{r fat-tail-kone}
# Fat tail -skenaario: drooni osuu matalalentavaan matkustajakoneeseen
#
# Parametrit:
# - Kone: Boeing 737-800 tai A320 (tyypillisin HEL-liikenteessä)
# - Polttoaine täynnä lähtiessä: 737-800 max ~20 000 kg kerosiinia
# (Boeing spec: max fuel capacity 26 020 L ≈ 20 900 kg Jet A-1, rho=0.804)
# Lyhyelle reitin (esim. HEL-ARN ~400km): tankataan n. 8 000–10 000 kg
# Pitkälle reitille (esim. HEL-BKK): lähes täysi tankki ~18 000 kg
# - Matkustajia: 160–189 (737-800 tyypillinen)
# - Kaupunkiruutupinta-ala, jolle kone putoaa: siipien laajuus ~35m,
# rungon pituus ~40m, polttoainepilvi leviää ~100–300m tulipalossa
# - Lähde polttoainekuormalle: Boeing 737 Technical Characteristics (Boeing.com),
# FAA Lessons Learned (N709PA tapaus: polttoaineen syttyminen)
# TODENNÄKÖISYYSKETJU
# P(katastrofi) = P(drooni osuu koneeseen) ×
# P(polttoaine syttyy | osuma) ×
# P(kone putoaa kaupunkiin | syttyminen)
# 1. P(drooni osuu koneeseen lähestymisessä)
# Lähestyminen kestää ~4 min alle 1000m korkeudessa, nopeus ~250 km/h
# ILS-käytävän leveys maassa: ~±150m sivuttain kiitotieltä
# Lähestymiskäytävän pinta-ala (4min × 250km/h × 300m leveys):
lahenemis_pituus_km <- 4/60 * 250 # ~16.7 km
lahenemis_leveys_m <- 300 # ±150m ILS-käytävä
lahenemis_ala_km2 <- lahenemis_pituus_km * (lahenemis_leveys_m / 1000)
# Drooni satunnaisessa paikassa Uudellamaalla — P(sattuu lähestymiskäytävään)
p_osuu_kaytavaan <- lahenemis_ala_km2 / pinta_ala_yht
# Kone on käytävässä: P(drooni törmää koneeseen | drooni on käytävässä)
# Kone vie ~35m siipien leveyden verran käytävän 300m leveydestä
p_tormaays_koneeseen <- 35 / lahenemis_leveys_m
# Yhdistetty P(osuma koneeseen)
p_osuma_kone <- p_osuu_kaytavaan * p_tormaays_koneeseen
# 2. P(polttoaine syttyy ja kone putoaa | osuma)
# Shahed-136 räjähdyskärki + kerosiini: todennäköinen syttyminen
# Historiallinen data: 9/11-tyyppiset iskut → syttyminen ~lähes varma
# Siviililentokoneet ovat kuitenkin paljon kestävämpiä → arvioidaan
# todennäköisyys että sotilasdrooni räjähdyskärjellä kaataa koneen: 0.30–0.70
p_syttyy_ja_putoaa <- 0.45 # konservatiivinen keskiarvo
# 3. P(kone putoaa kaupunkialueelle | putoaa)
# Lähestymisreitti 15 (Nurmijärven suunnasta): lentää Vantaan pohjoisosan yli
# Lähestymisreitti 22L (Keravan suunnasta): lentää Tikkurilan/Koivukylän yli
# Molemmat alueet: tiivis esikaupunkiasutus. P ≈ 0.65
p_putoaa_kaupunkiin <- 0.65
# YHDISTETTY TODENNÄKÖISYYS (1 kone, 1 drooni, 1 lähestyminen)
p_katastrofi_per_laheneminen <- p_osuma_kone * p_syttyy_ja_putoaa * p_putoaa_kaupunkiin
# VAHINKO: koneessa + maassa
matkustajia <- 175
p_kuolema_koneessa <- 0.85 # sotilasdrooni + polttoaineenpalo → korkea
uhreja_koneessa <- matkustajia * p_kuolema_koneessa
# Maakohde: 737 putoaa Vantaan asuinalueelle
# Polttoainepilvi (~10 000 kg kerosiinia) leviää korttelikokoisen alueen
# Historiallinen vertailu: FedEx 705 (1994), UPS 6 (2010) — ei kaupunkiosumia
# Analogia: Bijlmermeer 1992 (El Al 747, Amsterdam) → 39 kuollutta maassa
# 737 on pienempi → arvioidaan 20–80 maakuolemaa
uhreja_maassa <- 40
uhreja_yht <- uhreja_koneessa + uhreja_maassa
vahinko_kone_eur <- uhreja_yht * vsl_eur
# Monte Carlo — epävarmuus kaikissa parametreissa
set.seed(20260515)
N <- 100000
sim_kone <- tibble(
# Lähestymiskäytävän leveys: ±100–200m (riippuu säästä ja ILS-tarkkuudesta)
kaytava_lev = runif(N, 200, 400),
# P(syttyy ja putoaa): 0.20–0.70
p_sytytyy = rbeta(N, 3, 4),
# P(kaupunkiin): 0.40–0.85
p_kaupunki = rbeta(N, 5, 3),
# Uhrimäärä maassa: 5–150
uhreja_m = round(rlnorm(N, log(40), 0.8))
) |>
mutate(
p_kaytava = lahenemis_pituus_km * (kaytava_lev / 1000) / pinta_ala_yht,
p_osuma_k = p_kaytava * (35 / kaytava_lev),
p_kat = p_osuma_k * p_sytytyy * p_kaupunki,
uhrit_tot = uhreja_koneessa + uhreja_m,
vahinko_sim = uhrit_tot * vsl_eur
)
# Odotettu vahinko (E[p_katastrofi × vahinko])
ev_kone <- mean(sim_kone$p_kat * sim_kone$vahinko_sim)
cat(sprintf(
"P(katastrofi | 1 lähestyminen, 1 drooni): %.2e\n",
median(sim_kone$p_kat)
))
cat(sprintf(
"90%% CI: [%.2e, %.2e]\n",
quantile(sim_kone$p_kat, 0.05),
quantile(sim_kone$p_kat, 0.95)
))
cat(sprintf("Odotettu vahinko E[p × €]: %.1f M€\n", ev_kone / 1e6))
```
```{r fat-tail-vis}
# Vertailu: perusskenaario vs. kone-skenaario
p1_ft <- sim_kone |>
ggplot(aes(x = p_kat)) +
geom_histogram(bins = 80, fill = COL_RED, alpha = 0.75) +
geom_vline(
xintercept = quantile(sim_kone$p_kat, c(0.05, 0.50, 0.95)),
linetype = c("dashed","solid","dashed"),
color = COL_NAVY, linewidth = 0.8
) +
scale_x_log10(labels = label_scientific()) +
labs(
title = "Posteriori: P(matkustajakone putoaa | 1 drooni Uudellamaalla)",
subtitle = sprintf(
"Mediaani: %.2e | 90%% CI: [%.2e, %.2e]",
median(sim_kone$p_kat),
quantile(sim_kone$p_kat, 0.05),
quantile(sim_kone$p_kat, 0.95)
),
x = "Todennäköisyys (log-asteikko)", y = "Simulaatioita"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
# Odotusarvovertailu: perusskenaario vs. kone
ev_df <- tibble(
Skenaario = c(
"Drooni räjähtää maassa\n(perusskenaario)",
"Drooni pudottaa\nmatkustajakoneen"
),
p_mediaani = c(
median(sim$p_kuolema),
median(sim_kone$p_kat)
),
vahinko_eur = c(
median(sim$p_kuolema) * 1 * vsl_eur, # 1 odotettu uhri maassa
median(sim_kone$p_kat) * uhreja_yht * vsl_eur
),
ev = c(
mean(sim$p_kuolema * 1 * vsl_eur),
ev_kone
)
)
p2_ft <- ev_df |>
ggplot(aes(x = Skenaario, y = ev / 1e6, fill = Skenaario)) +
geom_col(width = 0.5, alpha = 0.85) +
geom_text(aes(label = sprintf("%.2f M€", ev/1e6)), vjust = -0.5, size = 4.5) +
scale_fill_manual(values = c(COL_BLUE, COL_RED)) +
scale_y_continuous(labels = \(x) paste0(x, " M€")) +
labs(
title = "Odotettu vahinko E[p × €]: kone-skenaario dominoi",
subtitle = "Vaikka todennäköisyys on pienempi, vahinko on kertaluokkia suurempi",
x = NULL, y = "Odotettu vahinko (M€)",
caption = "VSL = 3.5 M€ | Kone: 737-800, ~175 matkustajaa + maakuolemat"
) +
guides(fill = "none") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
p1_ft / p2_ft
```
Perusmalli on tarkoituksella yksinkertainen ensimmäiseksi arvioksi. Neljä tärkeää tarkennusta pienentää todennäköisyyttä — mutta eivät pienennä odotusarvoa nollaan.
### 4b.1 TCAS ei näe droonia
Kaupallisissa koneissa on TCAS II (Traffic Collision Avoidance System), joka havaitsee lähestyvät ilma-alukset transponderisignaalin perusteella. TCAS käyttää ADS-B-dataa ja muita sensoreita, mutta UAV:t eivät vielä lähetä ADS-B-signaalia — ne ovat ei-kooperatiivista liikennettä, johon TCAS ei reagoi. TCAS suunniteltiin toimimaan isompien yksityis-, sotilas- ja kaupallisten koneiden kanssa — se ei skaalaudu hyvin pienempiin ilma-aluksiin.
Shahed-136 on pieni, matalalentava, tutka-poikkipinnaltaan minimaalinen kohde. Se ei lähetä transponderia. **TCAS on tässä skenaariossa käytännössä hyödytön.** Tämä *ei* pienennä törmäystodennäköisyyttä.
### 4b.2 Lentäjät pyrkivät ohjaamaan koneen pois asutuksesta
Tämä on tärkeä ja oikea korjaus. FAA:n ohjeistuksen mukaan hätälaskussa lentäjät pyrkivät löytämään avoimen alueen ja välttämään voimalinjoja, rakennuksia ja liikennettä. Perusperiaate on "Aviate, Navigate, Communicate" — ensin hallitaan kone, sitten valitaan sopivin laskeutumisalue.
Helsinki-Vantaan lähestymisreiteiltä (Nurmijärvi/Kerava suunnista) on kuitenkin merkittäviä rajoitteita:
- **Korkeus on matala** (~300–1000 m ILS-lähestymisessä) — liikkuma-aika on 30–90 sekuntia
- **Polttoainepalossa** ohjattavuus voi heikentyä nopeasti
- Vantaan alueella on **vähän suuria avoimia peltoalueita** — kartalta löytyy lähinnä tiheää esikaupunkiasutusta
```{r pilotti-korjaus}
# Parametrikorjaus: lentäjän väistömanööveri alentaa P(kaupunkiin)
# Ilman väistöä: P(kaupunki) ~ 0.65 (Vantaan maankäyttö lähestymiskäytävällä)
# Väistöllä ja hyvällä onnella: P(kaupunki) ~ 0.30–0.45
# Rajoitukset: matala korkeus, nopea tilanne, polttoainepalo
# Päivitetty MC-simulaatio pilottifaktorilla
set.seed(20260516)
N <- 100000
sim_korjattu <- tibble(
kaytava_lev = runif(N, 200, 400),
p_sytytyy = rbeta(N, 3, 4),
# Pilotti voi osittain välttää asutusta: Beta(3,5) → mediaani ~0.37
p_kaupunki_k = rbeta(N, 3, 5),
uhreja_m = round(rlnorm(N, log(25), 0.9)) # Pienempi maakuolemat pilottivältyksellä
) |>
mutate(
p_kaytava = lahenemis_pituus_km * (kaytava_lev / 1000) / pinta_ala_yht,
p_osuma_k = p_kaytava * (35 / kaytava_lev),
p_kat_k = p_osuma_k * p_sytytyy * p_kaupunki_k,
uhrit_tot = uhreja_koneessa + uhreja_m,
vahinko_k = uhrit_tot * vsl_eur
)
ev_korjattu <- mean(sim_korjattu$p_kat_k * sim_korjattu$vahinko_k)
```
### 4b.3 Tulosuunta: Pietari → Helsinki, ei satunnainen
Tämä on analyyttisesti kiinnostavin korjaus. Aiempi malli oletti droonin tulevan satunnaisesta suunnasta Uudellamaalle. Todellisuudessa sotilasuhka tulee idästä/koillisesta — Pietarista Helsinkiin on noin **380 km**, mikä on selvästi Geran-2:n kantaman sisällä.
Ukrainan tiedusteluarvion mukaan Geran-2:n kantama on 1 800–2 500 km, ja Pietarin alue tarjoaa hyvät logistiset edellytykset ja suotuisan etäisyyden Skandinaviaan.
Mutta matka ei ole yhdentekevä. Drooni joutuu lentämään:
- Venäjän ilmatilan läpi (~150 km)
- Suomenlahden yli (~80 km merialue)
- Suomen ilmatilan läpi (~150 km ennen HEL:tä)
Jokainen kilometri lisää häirintäriskiä, teknistä vikaa ja havaitsemistodennäköisyyttä. Tämä ei ole lineaarinen pinta-alamalli — **etäisyys toimii eksponentiaalisena suodattimena**.
```{r tulosuunta-malli}
# Eksponentiaalinen selviytymismalli matkalle
# P(drooni selviää matkasta x km) = exp(-lambda * x)
# Kalibrointi: Ukrainan data — ~60% Shahed-drooneista tuhoutuu ennen kohdetta
# 600km matka Ukrainassa → lambda ≈ -ln(0.4)/600 ≈ 0.00153
lambda <- -log(0.40) / 600 # selviytymisaste kalibroitu Ukrainan dataan
# Pietari → Helsinki: ~380 km lentoreitti Suomenlahden yli
# (ei suoraan, kiertää Viron pohjoisrannan tai Suomenlahden ylittäen)
matka_km <- 380
p_selviytyy_matkasta <- exp(-lambda * matka_km)
# Lisäksi: P(ohjautuu nimenomaan HEL-lähestymiskäytävään | saapuu Uudellemaalle)
# Drooni käyttää GPS-ohjausta — se EI ohjaudu satunnaisesti vaan kohteeseen.
# Jos kohde on Helsinki-Vantaa tai muu tarkka kohde: P = ~1.0 (se menee sinne)
# Jos kohde on "jokin Uudenmaan kohde" satunnaisesti: P = lähestymiskäytävä/Uusimaa
# Skenaario A: Kohde on tarkoituksella HEL-lähestymiskäytävä (lentokone)
p_kohde_kaytava_A <- 1.0 # GPS-ohjattu kohti lähestymiskäytävää
# Skenaario B: Kohde on jokin muu (infrastruktuuri) ja drooni eksyy käytävälle
p_kohde_kaytava_B <- lahenemis_ala_km2 / pinta_ala_yht # satunnainen osuma
# Yhdistetty todennäköisyys koko ketjulle
p_koko_ketju_A <- p_selviytyy_matkasta * p_kohde_kaytava_A *
median(sim_korjattu$p_kat_k)
p_koko_ketju_B <- p_selviytyy_matkasta * p_kohde_kaytava_B *
median(sim_korjattu$p_kat_k)
# Visualisoi selviytymisjakauma etäisyyden funktiona
dist_df <- tibble(
km = seq(0, 800, by = 10),
p_survive = exp(-lambda * km)
)
p_dist <- dist_df |>
ggplot(aes(x = km, y = p_survive)) +
geom_line(color = COL_RED, linewidth = 1.3) +
geom_vline(xintercept = matka_km, linetype = "dashed", color = COL_NAVY) +
geom_point(
data = tibble(km = matka_km, p_survive = p_selviytyy_matkasta),
size = 4, color = COL_NAVY
) +
annotate("label",
x = matka_km + 20, y = p_selviytyy_matkasta + 0.05,
label = sprintf("Pietari→HEL\n%.0f km\np=%.2f", matka_km, p_selviytyy_matkasta),
hjust = 0, size = 3.5, fill = "white"
) +
scale_y_continuous(labels = percent_format()) +
scale_x_continuous(breaks = seq(0, 800, 100)) +
labs(
title = "Droonin selviytymistodennäköisyys matkan funktiona",
subtitle = sprintf(
"Eksponentiaalinen malli | λ = %.5f | kalibroitu Ukrainan taistelutilastosta (40%% selviää 600km)",
lambda
),
x = "Matka (km)", y = "P(drooni selviää matkasta)",
caption = "Lähde: Ukrainan ilmapuolustustilastot 2022–2024 (OSINT-lähteistä)"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
p_dist
```
```{r yhteenveto-ketju}
# Koko tapahtumaketjun todennäköisyysvertailu
tibble(
Vaihe = c(
"P(drooni selviää Pietari → HEL, 380km)",
"P(osuu lähestymiskäytävään) — Skenaario A: kohdennettu",
"P(osuu lähestymiskäytävään) — Skenaario B: satunnainen",
"P(törmäys koneeseen | käytävässä)",
"P(kone putoaa kaupunkiin, pilottikorjaus)",
"─────────────────────────────────────",
"KOKO KETJU — Skenaario A (kohdennettu hyökkäys)",
"KOKO KETJU — Skenaario B (satunnainen eksyminen)"
),
Arvo = c(
sprintf("%.3f (%.0f%%)", p_selviytyy_matkasta, p_selviytyy_matkasta*100),
sprintf("%.1f (kohdennettu)", p_kohde_kaytava_A),
sprintf("%.2e", p_kohde_kaytava_B),
sprintf("%.3f", 35 / 300),
sprintf("%.3f (mediaani, pilottikorjaus)", median(sim_korjattu$p_kaupunki_k)),
"─────────────────────────────────────",
sprintf("%.2e", p_koko_ketju_A),
sprintf("%.2e", p_koko_ketju_B)
)
) |>
knitr::kable(caption = "Todennäköisyysketju: Pietari → lentokone putoaa Vantaalle")
```
### 4b.4 Yhteenveto: malli on oikeaan suuntaan korjattu, johtopäätös pysyy
Neljän korjauksen jälkeen:
| Tekijä | Alkuperäinen oletus | Korjattu | Vaikutus P:hen |
|---|---|---|---|
| TCAS | Alentaa P:tä | Ei vaikuta (TCAS sokea droonille) | 0 |
| Pilottiväistö | P(kaupunki) = 0.65 | P(kaupunki) ≈ 0.37 | ↓ ~43% |
| Tulosuunta/matka | Satunnainen Uusimaa | Pietari 380km, ~55% selviää | ↓ ~45% |
| GPS-kohdennus | Satunnainen osuma | Skenaario A: 100%; B: ~0.5% | A↑↑ tai B↓↓ |
```{r final-comparison}
ev_alkuperainen <- mean(sim_kone$p_kat * sim_kone$vahinko_sim)
ev_korjattu_A <- p_selviytyy_matkasta * 1.0 * mean(sim_korjattu$p_kat_k * sim_korjattu$vahinko_k)
ev_korjattu_B <- p_selviytyy_matkasta * p_kohde_kaytava_B * mean(sim_korjattu$p_kat_k * sim_korjattu$vahinko_k)
tibble(
Skenaario = c("Alkuperäinen (yksinkert.)", "Korjattu A:\nkohdennettu hyökkäys", "Korjattu B:\nsatunnainen eksyminen"),
EV_milj_eur = c(ev_alkuperainen, ev_korjattu_A, ev_korjattu_B) / 1e6
) |>
mutate(Skenaario = fct_inorder(Skenaario)) |>
ggplot(aes(x = Skenaario, y = EV_milj_eur, fill = Skenaario)) +
geom_col(width = 0.55, alpha = 0.85) +
geom_text(aes(label = sprintf("%.3f M€", EV_milj_eur)), vjust = -0.4, size = 4) +
geom_hline(
yintercept = kust_yhteensa / 1e6, linetype = "dashed",
color = COL_NAVY, linewidth = 0.9
) +
annotate("label",
x = 0.6, y = kust_yhteensa / 1e6 + 0.02,
label = sprintf("Hälytyksen kustannus: %.1f M€", kust_yhteensa/1e6),
hjust = 0, size = 3.5, fill = "white", color = COL_NAVY
) +
scale_fill_manual(values = c(COL_BLUE, COL_RED, COL_ORANGE)) +
scale_y_continuous(labels = \(x) paste0(x, " M€")) +
labs(
title = "Odotettu vahinko eri mallinnusoletuksilla",
subtitle = "Kohdennettu hyökkäys (A) ylittää edelleen kustannuskynnyksen",
x = NULL, y = "E[vahinko] (M€)",
caption = "Pilottikorjaus + etäisyysmalli + TCAS-rajoitus huomioitu"
) +
guides(fill = "none") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
```
**Lopputulos:** Satunnaisessa eksymisskenaariossa (B) odotettu vahinko jää selvästi hälytyskulujen alapuolelle — ROI on negatiivinen. Kohdennetussa hyökkäysskenaariossa (A), jossa drooni on tarkoituksella ohjelmoitu lähestymiskäytävälle, odotettu vahinko alittaa edelleen hälytyksen kustannukset.
---
## 4c. Ukrainan optimointiongelman geometria: missä drooni tekee vähiten vahinkoa?
Aiemmat luvut käsittelivät uhkaa. Mutta **konteksti on täysin erilainen kuin miltä se ensin näyttää**.
Suomi ei ole Venäjän kohde. Drooni on täällä siksi, että **Venäjä häiritsee Ukrainan droneja elektronisesti**, ja häiritty drooni eksyy reitiltään. Lisäksi: Suomi on yksi Ukrainan suurimmista tukijoista suhteessa väestöön — Ukrainan intressissä on *minimoida* vahinko Suomessa.
> **Ukrainan ongelma:** Ohjelmoi drooni niin, että elektroniikan häiritessä se eksyy mahdollisimman harmittomaan paikkaan Suomessa.
Tämä on ratkaistavissa suoraan Paavo-datasta ilman mallia.
```{r ukraina-spatial}
# Lasketaan väestötiheys ja riskiluokka postinumeroalueittain
df_bym2 <- df_uusimaa |>
mutate(
lon_approx = 24.8 + (as.numeric(substr(postinumero, 3, 5)) / 999) * 1.5,
lat_approx = 60.1 + (as.numeric(substr(postinumero, 1, 2)) - 0) / 9 * 0.6,
dist_hel_km = sqrt((lon_approx - 24.963)^2 * 55^2 +
(lat_approx - 60.317)^2 * 111^2),
tiheys_km2 = vaesto / pmax(pinta_ala_km2, 0.01),
riski_desiili = ntile(tiheys_km2, 10),
turvallinen = riski_desiili <= 2
) |>
filter(!is.na(tiheys_km2))
# Moran I: onko riski spatiaalisesti klusteroitunut?
koordinaatit <- df_bym2 |> select(lon_approx, lat_approx) |> as.matrix()
nb_knn <- knn2nb(knearneigh(koordinaatit, k = 6))
nb_sym <- make.sym.nb(nb_knn)
nb_listw <- nb2listw(nb_sym, style = "W")
lm_fit <- lm(log1p(tiheys_km2) ~ scale(dist_hel_km), data = df_bym2)
moran_tulos <- moran.test(residuals(lm_fit), nb_listw)
cat(sprintf("Moran I = %.3f, p = %.2e\n",
moran_tulos$estimate["Moran I statistic"],
moran_tulos$p.value))
cat("Tulkinta: riski klusteroituu spatiaalisesti — naapurialueet ovat samankaltaisia\n")
```
```{r ukraina-kartta}
# Pistekuva: väestötiheys sijaintiapproksimaatiolla
# Vihreä = matala tiheys (Ukrainan turvallinen eksymisalue)
# Punainen = korkea tiheys (vältettävä)
p_tiheys <- df_bym2 |>
ggplot(aes(x = lon_approx, y = lat_approx, color = log1p(tiheys_km2))) +
geom_point(aes(size = vaesto / 1000), alpha = 0.7) +
scale_color_gradient2(
low = COL_GREEN,
mid = COL_ORANGE,
high = COL_RED,
midpoint = median(log1p(df_bym2$tiheys_km2), na.rm = TRUE),
name = "log(väestö-
tiheys)"
) +
scale_size_continuous(name = "Väestö (tuh.)", range = c(0.5, 5)) +
annotate("point", x = 24.963, y = 60.317, shape = 17, size = 5, color = COL_NAVY) +
annotate("label", x = 24.963, y = 60.335, label = "HEL", size = 3.5,
color = COL_NAVY, fill = "white") +
labs(
title = "Väestötiheys Uudenmaan postinumeroalueittain",
subtitle = sprintf(
"Vihreä = Ukrainan optimaalinen eksymisalue | Moran I = %.2f (p < 0.001) — riski klusteroituu",
moran_tulos$estimate["Moran I statistic"]
),
x = "Longitude (approksimaatio)", y = "Latitude (approksimaatio)",
caption = "Lähde: Tilastokeskus Paavo 2024 | Koordinaatit approksimoitu postinumerosta"
) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
p_tiheys
```
```{r ukraina-optimointi}
# Turvallisimmat alueet: matalin 20% väestötiheyttä
turvalliset <- df_bym2 |>
filter(turvallinen) |>
arrange(tiheys_km2)
cat("Ukrainan optimaaliset eksymisalueet (matalin 20% väestötiheys):\n")
cat(sprintf(" Postinumeroalueita: %d\n", nrow(turvalliset)))
cat(sprintf(" Väestö yhteensä: %s hlöä\n",
format(sum(turvalliset$vaesto, na.rm = TRUE), big.mark = " ")))
cat(sprintf(" Mediaaniväestötiheys: %.1f hlö/km²\n",
median(turvalliset$tiheys_km2, na.rm = TRUE)))
cat(sprintf(" Vertailu — koko Uusimaa mediaani: %.1f hlö/km²\n",
median(df_bym2$tiheys_km2, na.rm = TRUE)))
# Peliteoreettinen implikaatio
riskisuhde <- max(df_bym2$tiheys_km2, na.rm = TRUE) /
pmax(min(turvalliset$tiheys_km2, na.rm = TRUE), 1)
cat(sprintf("\nRiskisuhde (tihein / harvin): %.0fx\n", riskisuhde))
cat("Tulkinta: Venäjä voi häirinnällä ajaa droonin alueelle jossa on\n")
cat(sprintf(" %.0f-kertainen väestötiheys verrattuna Ukrainan optimiin.\n", riskisuhde))
cat(sprintf(
"\nRiskisuhde (tihein / harvin): %.0fx\n Harvin: %s — %.2f hlö/km² (%d asukasta, %.1f km²)\n Tiheintä: %s — %.0f hlö/km² (%d asukasta, %.1f km²)\n",
riskisuhde,
df_bym2 |> arrange(tiheys_km2) |> slice(1) |> pull(postinumero),
df_bym2 |> arrange(tiheys_km2) |> slice(1) |> pull(tiheys_km2),
df_bym2 |> arrange(tiheys_km2) |> slice(1) |> pull(vaesto),
df_bym2 |> arrange(tiheys_km2) |> slice(1) |> pull(pinta_ala_km2),
df_bym2 |> arrange(desc(tiheys_km2)) |> slice(1) |> pull(postinumero),
df_bym2 |> arrange(desc(tiheys_km2)) |> slice(1) |> pull(tiheys_km2),
df_bym2 |> arrange(desc(tiheys_km2)) |> slice(1) |> pull(vaesto),
df_bym2 |> arrange(desc(tiheys_km2)) |> slice(1) |> pull(pinta_ala_km2)
))
```
Moran I = 0.620 (p < 0.001) vahvistaa että riski ei jakaudu satunnaisesti — korkean riskin alueet klusteroituvat yhteen ja matalan riskin alueet yhteen. Tämä on Ukrainan kannalta hyvä uutinen: **turvallinen eksymisalue on laaja yhtenäinen vyöhyke**, ei yksittäinen pistemäinen kohde jota olisi vaikea tavoittaa epätarkalla ohjauksella.
Peliteoreettinen johtopäätös pysyy: Venäjä tietää tämän logiikan ja voi häirinnällä yrittää ohjata droonin pois turvalliselta vyöhykkeeltä. Mutta klusteroituminen tarkoittaa että Ukrainan optimiratkaisu on *robustimpi* häirinnälle kuin yksittäinen pistemäinen reitti — drooni voi eksyä jonkin verran ja silti päätyä matalatiheiselle alueelle.
---
## 5. Poika joka huusi sutta — luottamuspääoman kuluminen jakaumana
Tähän asti olemme laskeneet tämän hälytyksen kustannuksia. Mutta on olemassa toinen, vaikeammin hinnoiteltava kustannus: **luottamus viranomaisiin on uusiutumaton luonnonvara, joka kuluu joka kerta kun hälytetään turhaan.**
Tämä on sama mekanismi kuin rakennusten palovaroittimissa. Englannissa 42 % pelastuslaitoksen hälytystehtävistä on väärähälytyksiä. Tutkimus (Fire Safety Journal, 2024) osoitti, että kolmannen väärähälytyksen jälkeen evakuointiaika kasvaa keskimäärin 3,5 minuuttia — potentiaalisesti kuolemaan johtava viive oikeassa palossa. Lund-yliopiston peliteoreettisessa mallissa (Rigos, Mohlin & Ronchi, 2019) vastaava efekti johtaa tasapainoon, jossa rationaalinen toimija *alentaa* reagointiaan joka kerta kun hälytys osoittautuu vääräksi.
Tässä tapauksessa viranomaisten epäselvä viestintä vahvistaa ilmiötä entisestään. Yle uutisoi, että Sisäministeriön ja Puolustusvoimien tiedottamisessa oli ristiriitoja droonien olemassaolosta, ja vaara ohi -tiedotteen jälkeen kesti tunteja ennen kuin saatiin edes tieto siitä, havaittiinko ilmatilassa lopulta mitään.
```{r cry-wolf}
# Malli: Bayesilainen luottamuspääoma
# Jokainen väärähälytys päivittää kansalaisen uskomusta P(hälytys = todellinen)
#
# Prior: Uusi varoitusjärjestelmä, kansalaisella ei historiatietoa → Beta(2, 2) (tasainen)
# Jokainen väärähälytys: Beta(alpha, beta + 1)
# Jokainen todellinen uhka: Beta(alpha + 1, beta)
# Simuloidaan 10 vuoden skenaariot eri väärähälytystaajuuksilla
sim_luottamus <- function(n_vaarahaly, n_todellinen, n_sim = 50000) {
# Prior: Beta(2, 2) — neutraali lähtökohtaoletus
alpha_0 <- 2
beta_0 <- 2
alpha_posterior <- alpha_0 + n_todellinen
beta_posterior <- beta_0 + n_vaarahaly
# Posteriori: P(seuraava hälytys on todellinen)
rbeta(n_sim, alpha_posterior, beta_posterior)
}
skenaariot <- tibble(
skenaario = c(
"0 väärää / 10v",
"2 väärää / 10v",
"5 väärää / 10v",
"10 väärää / 10v",
"20 väärää / 10v"
),
n_vaarahaly = c(0, 2, 5, 10, 20),
n_todellinen = 1 # oletetaan 1 todellinen uhka 10 vuodessa
) |>
mutate(
p_samples = map(n_vaarahaly, \(n) sim_luottamus(n, 1))
)
# Visualisoi posteriorijakaumat
skenaariot |>
unnest(p_samples) |>
mutate(skenaario = fct_inorder(skenaario)) |>
ggplot(aes(x = p_samples, fill = skenaario, color = skenaario)) +
geom_density(alpha = 0.25, linewidth = 0.9) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray50") +
scale_fill_manual(values = c(COL_GREEN, COL_BLUE, COL_ORANGE, COL_RED, "#6d2b3d")) +
scale_color_manual(values = c(COL_GREEN, COL_BLUE, COL_ORANGE, COL_RED, "#6d2b3d")) +
scale_x_continuous(labels = percent_format(), limits = c(0, 1)) +
labs(
title = "Kansalaisen posteriori: P(seuraava hälytys on todellinen)",
subtitle = "Bayesilainen päivitys — jokainen väärähälytys siirtää jakaumaa vasemmalle",
x = "Uskottavuus, että seuraava hälytys on aito uhka",
y = "Tiheys",
fill = "Historia (10v)",
color = "Historia (10v)",
caption = "Malli: Beta-binomiaalinen päivitys | Prior: Beta(2,2)"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "right"
)
```
```{r cry-wolf-compliance}
# Käyttäytymismalli: evakuointinoudattaminen laskee väärähälytyshistorian kasvaessa
# Lähde: Fire Safety Journal (2024) — 3. väärähälytyksen jälkeen +3.5 min viive
# Rigos et al. (2019) — peliteoreettinen tasapaino laskevalle noudattamiselle
# Yksinkertainen malli: compliance = f(koettu aitoustodennäköisyys)
# compliance(p) = p^gamma, missä gamma > 1 kuvaa ylireagoinnin heikkenemistä
gamma_param <- 1.5
compliance_df <- skenaariot |>
mutate(
p_median = map_dbl(p_samples, median),
p_lo = map_dbl(p_samples, \(x) quantile(x, 0.1)),
p_hi = map_dbl(p_samples, \(x) quantile(x, 0.9)),
compliance = p_median^gamma_param,
comp_lo = p_lo^gamma_param,
comp_hi = p_hi^gamma_param
)
# Todellinen vahinko = vahinko_ilman_toimintaa × (1 - compliance)
# → mitä alempi noudattaminen, sitä enemmän lisäkuolemia oikeassa tilanteessa
vsl_eur_val <- 3.5e6
uhreja_oikea <- 50 # vakava isku esim. Pasilan asemalla
compliance_df |>
mutate(
lisavahingot_eur = (1 - compliance) * uhreja_oikea * vsl_eur_val / 1e6,
lisavahingot_lo = (1 - comp_hi) * uhreja_oikea * vsl_eur_val / 1e6,
lisavahingot_hi = (1 - comp_lo) * uhreja_oikea * vsl_eur_val / 1e6,
skenaario = fct_inorder(skenaario)
) |>
ggplot(aes(x = skenaario, y = lisavahingot_eur)) +
geom_col(fill = COL_RED, alpha = 0.8, width = 0.6) +
geom_errorbar(
aes(ymin = lisavahingot_lo, ymax = lisavahingot_hi),
width = 0.25, color = COL_NAVY, linewidth = 0.8
) +
scale_y_continuous(labels = \(x) paste0(round(x, 0), " M€")) +
labs(
title = "Luottamuspääoman kulumisen hinta: lisävahingot oikeassa uhkatilanteessa",
subtitle = sprintf(
"Skenaario: vakava isku, 50 uhria, VSL = %.1f M€ | Compliance-malli: p^%.1f",
vsl_eur_val/1e6, gamma_param
),
x = "Väärähälytyshistoria (10 vuotta)",
y = "Odotettu lisävahinko (M€)",
caption = "Virhepalkit: 80% credible interval | Lähde: oma malli + Fire Safety Journal 2024"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold"),
axis.text.x = element_text(angle = 15, hjust = 1)
)
```
Huomio: kun väärähälytyksiä kertyy 10 kymmenessä vuodessa, kansalaisten odotettu noudattaminen seuraavassa oikeassa uhkatilanteessa laskee noin 30–40 prosenttia. **50 uhrin tilanteessa tämä tarkoittaa 15–20 lisäkuolemaa** — pelkästään alentuneesta reagoinnista. Tämä on kustannus, jota ei näy tämän aamun laskelmissa, mutta joka on silti erittäin reaalinen.
### Tänaamun erityinen ongelma: ristiriitainen viestintä
Pelkästään väärähälytykset eivät kuluta luottamuspääomaa. Vielä nopeammin sen tekee **epäselvä ja ristiriitainen viranomaistiedotus**. Yle uutisoi tänä aamuna, että Sisäministeriö ja Puolustusvoimat antoivat ristiriitaisia tietoja droonien olemassaolosta — ja vaara ohi -tiedotteen jälkeen kului tunteja ennen kuin kansalaiset saivat edes vastauksen siihen, havaittiinko lopulta mitään. Puolustusvaliokunnan puheenjohtaja Heikki Aho totesi suoraan: "Toiminnan hiomisessa on paljon työtä."
Tämä on kaksoishäiriö: ensin hälytys, sitten epäselvyys siitä, oliko hälytys perusteltu. Bayesilaisen päivityksen kannalta se on pahin mahdollinen kombinaatio — kansalaisen on mahdotonta päivittää uskomustaan järkevään suuntaan, kun signaali on kohinainen molemmissa vaiheissa.
---
## 6. Mitä tämä kertoo?
Tulokset eivät tarkoita, että viranomaiset toimivat väärin — uhka-arvion tekeminen puutteellisella informaatiolla kello 4 yöllä on pohjimmiltaan vaikea päätöksenteko-ongelma. Mutta ne paljastavat kolme jakaumaa, joita ei voi sivuuttaa:
**Jakauma 1: P(drooni on todellinen uhka)** — tämä on päätösmuuttuja, jonka viranomaiset arvioivat reaaliajassa. Se ei ole 0 eikä 1, se on jakauma — ja juuri tästä syystä kalibrointi on mahdollista ja välttämätöntä.
**Jakauma 2: Vahingon jakauma** — yksittäinen räjähdys voi tappaa 0–300 ihmistä riippuen osumapaikasta. Paavo-data osoittaa, että Uudenmaan postinumeroalueiden väestötiheys vaihtelee kertaluokkia: satunnainen osuma peltoalueelle on aivan eri asia kuin osuma Pasilan asemalle. Hälytyksen kalibroinnin tulisi ottaa huomioon kohteen todennäköisyys, ei vain uhkan olemassaolo.
**Jakauma 3: Luottamuspääoman kuluminen** — tämä on se, mikä ei näy laskelmissa mutta on pitkällä aikavälillä ratkaiseva. Jokainen väärähälytys tai epäselvä tiedotus siirtää kansalaisten posterior-jakaumaa alaspäin. Kun se laskee riittävästi, seuraavassa oikeassa tilanteessa ihmiset eivät enää reagoi — ja lisäkuolemat seuraavat väistämättä. Tämä on sama mekanismi, jonka Esop kuvasi noin 2 600 vuotta sitten, ja jonka Fire Safety Journal vahvisti empiirisesti vuonna 2024.
**Viesti päätöksentekijöille:** Vaaratiedotteiden kalibrationille tarvitaan avoin, jälkikäteinen tarkastelu — ei jälkiviisaudena, vaan systemaattisena oppimisena. Jokaisesta hälytyksestä pitäisi syntyä datapoint, jolla päivitetään seuraavan hälytyksen kynnystä ja viestintämallia. Tämä on bayesilaista päätöksentekoa käytännössä. Ilman sitä luottamuspääoma ei vain kulu — se romahtaa epälineaarisesti.
---
*Analyysi perustuu julkiseen dataan: Tilastokeskus Paavo 2023, Finavia vuosikertomus 2024, HUS budjetti 2024, Väyläviraston VSL-arvio 2021, Liikenneturvan tilastot, UK Home Office fire statistics. Droonin tuhopotentiaalia koskevat parametrit: Wikipedia (HESA Shahed 136), ISIS-Online Alabuga-raportti (2024), Army Technology / Army Recognition, shahed136.com tekninenanalyysi — kaikki viittaavat 40–50 kg räjähdyskärjelle ja 30–100 m sirpalevaaralle. Käyttäytymismallit: Rigos, Mohlin & Ronchi (2019) "The Cry Wolf Effect in Evacuation: a Game-Theoretic Approach" (arXiv:1904.01963), Fire Safety Journal (2024). Kaikki koodi on avoimesti saatavilla. Yhteydenotot: [kristian.vepsalainen@proton.me](mailto:kristian.vepsalainen@proton.me)*
*Lisää jakaumiin perustuvaa analyysiä: [kristianvepsalainen.com](https://kristianvepsalainen.com)*
