预测核辐射强度的分布
空间数据分析
黄湘云
统计之都 @ 北京会议中心
2023年7月12日
目录
核污染研究的背景和意义
数据探索和分析的过程
- 描述数据的空间分布
- 广义线性模型及模型诊断
- 空间广义线性混合效应模型
- 模型拟合及准备预测位置
- 模型预测及展示预测分布
氢弹核试验
这颗氢弹的当量相当于广岛原子弹的 1000 倍。 由于核辐射扩散，导致朗格拉普环礁的主岛受到核辐射1。
朗格拉普岛
图 1: 朗格拉普岛地理位置特殊，是一个重要的军事基地，建有环岛公路和机场。卫星拍摄的影像数据（2020 年）来自 https://rmi-data.sprep.org/
采集的核辐射数据（部分）
表格 1: 采集的核辐射数据
| 横坐标 | 纵坐标 | 粒子数目 | 统计时间 |
|---|---|---|---|
| -6050 | -3270 | 75 | 300 |
| -6050 | -3165 | 371 | 300 |
| -5925 | -3320 | 1931 | 300 |
| -5925 | -3165 | 4357 | 300 |
| -5800 | -3350 | 2114 | 300 |
| -5800 | -3165 | 2318 | 300 |
| 横坐标 | 纵坐标 |
|---|---|
| -5509.236 | -3577.438 |
| -5544.821 | -3582.250 |
| -5561.604 | -3576.926 |
| -5580.780 | -3574.535 |
| -5599.687 | -3564.288 |
| -5605.922 | -3560.910 |
采样数据的空间分布
岛上各采样点的核辐射强度
图 4: 岛上各采样点的核辐射强度
广义线性模型及拟合
核辐射是由放射元素衰变产生的，通常用单位时间释放出来的粒子数目表示辐射强度，因此，建立如下泊松型广义线性模型来拟合核辐射强度。
log(λi)yi=β∼Poisson(tiλi)
待估参数 β 。
fit_rongelap_poisson <- glm(counts ~ 1,
family = poisson(link = "log"),
offset = log(time), data = rongelap
)
summary(fit_rongelap_poisson)
Call:
glm(formula = counts ~ 1, family = poisson(link = "log"), data = rongelap,
offset = log(time))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.013954 0.001454 1385 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 61567 on 156 degrees of freedom
Residual deviance: 61567 on 156 degrees of freedom
AIC: 63089
Number of Fisher Scoring iterations: 4模型诊断 I
图 5: 残差的空间分布
颜色深和颜色浅的点分别聚集在一起，且与周围点的颜色呈现层次变化，拟合残差存在一定的空间相关性。
模型诊断 II
空间广义线性混合效应模型
广义线性模型并没有考虑距离相关性，它认为各个观测点的数据是相互独立的。因此，在广义线性模型的基础上添加位置相关的随机效应，用以刻画未能直接观测到的潜在影响。根据 137Cs 放出伽马射线，在 n=157 个采样点，分别以时间间隔 ti 测量辐射量 y(xi)，建立泊松型空间广义线性混合效应模型(Diggle, Tawn, 和 Moyeed 1998)。
log{λ(xi)}y(xi)=β+S(xi)∼Poisson(tiλ(xi))(1)
其中，β 表示截距，相当于平均水平，λ(xi) 表示位置 xi 处的辐射强度，S(xi) 表示位置 xi 处的空间效应，S(x),x∈D⊂R2 是二维平稳空间高斯过程 S 的具体实现。 D 表示研究区域，可以理解为朗格拉普岛，它是二维实平面 R2 的子集。
随机过程的自协方差函数
随机过程 S(x) 的自协方差函数常用的有指数型、幂二次指数型（高斯型）和梅隆型，形式如下：
Cov{S(xi),S(xj)}Cov{S(xi),S(xj)}Cov{S(xi),S(xj)}Kν(x)=σ2exp(−ϕ∥xi−xj∥2)=σ2exp(−2ϕ2∥xi−xj∥22)=σ2Γ(ν)21−ν(2νϕ∥xi−xj∥2)νKν(2νϕ∥xi−xj∥2)=∫0∞exp(−xcosht)cosh(νt)dt
待估参数：代表方差的 σ2 和代表范围的 ϕ 。当 ν=1/2 时，梅隆型退化为指数型。
模型拟合
library(spaMM)
# 对数高斯分布作为参考
fit_rongelap_gaussian <- fitme(
log(counts / time) ~ 1 + Matern(1 | cX + cY),
data = rongelap, fixed = list(nu = 0.5), method = "REML"
)
# 泊松分布
fit_rongelap_poisson <- fitme(
formula = counts ~ 1 + Matern(1 | cX + cY) + offset(log(time)),
family = poisson(link = "log"), data = rongelap,
fixed = list(nu = 0.5), method = "REML"
)
summary(fit_rongelap_poisson)formula: counts ~ 1 + Matern(1 | cX + cY) + offset(log(time))
Estimation of corrPars and lambda by REML (p_bv approximation of restricted logL).
Estimation of fixed effects by ML (p_v approximation of logL).
Estimation of lambda by 'outer' REML, maximizing restricted logL.
family: poisson( link = log )
------------ Fixed effects (beta) ------------
Estimate Cond. SE t-value
(Intercept) 1.829 0.08797 20.78
--------------- Random effects ---------------
Family: gaussian( link = identity )
--- Correlation parameters:
1.nu 1.rho
0.500000000 0.009211777
--- Variance parameters ('lambda'):
lambda = var(u) for u ~ Gaussian;
cX + cY : 0.3069
# of obs: 157; # of groups: cX + cY, 157
------------- Likelihood values -------------
logLik
logL (p_v(h)): -1318.010
Re.logL (p_b,v(h)): -1319.522采用 spaMM 包(Rousset 和 Ferdy 2014) 拟合模型。输出结果中，固定效应中的 (Intercept) 对应模型参数 β=1.829 ，随机效应中的方差参数 lambda 对应模型参数 σ2=0.3069 ，随机效应中相关参数 1.rho 的倒数对应模型范围参数 ϕ=1/0.00921=108.58 。
准备数据：构造网格
代码
library(sf)
rongelap_sf <- st_as_sf(rongelap, coords = c("cX", "cY"), dim = "XY")
rongelap_coastline_sf <- st_as_sf(rongelap_coastline, coords = c("cX", "cY"), dim = "XY")
# 点转线
rongelap_coastline_sfp <- st_cast(st_combine(st_geometry(rongelap_coastline_sf)), "POLYGON")
# 缓冲区
rongelap_coastline_buffer <- st_buffer(rongelap_coastline_sfp, dist = 50)
# 矩形网格
rongelap_coastline_grid <- st_make_grid(rongelap_coastline_buffer, n = c(150, 75))
图 9: 网格操作
采用 sf 包(Pebesma 2018)和 stars 包(Pebesma 和 Bivand 2023)操作空间数据。
准备数据：未采样的位置
代码
# 将 sfc 类型转化为 sf 类型，准备取交集
rongelap_coastline_grid <- st_as_sf(rongelap_coastline_grid)
rongelap_coastline_buffer <- st_as_sf(rongelap_coastline_buffer)
# 边界约束内的网格
rongelap_grid <- rongelap_coastline_grid[rongelap_coastline_buffer, op = st_intersects]
# 计算网格中心点坐标
rongelap_grid_centroid <- st_centroid(rongelap_grid)
图 10: 朗格拉普岛规则网格划分结果
模型预测
# sf 类型转 data.frame 类型
rongelap_grid_df <- as.data.frame(st_coordinates(rongelap_grid_centroid))
colnames(rongelap_grid_df) <- c("cX", "cY")
rongelap_grid_df$time <- 1
# 对数高斯分布
# 预测
rongelap_grid_pred <- predict(fit_rongelap_gaussian,
newdata = rongelap_grid_df, type = "response"
)
rongelap_grid_df$pred_sk <- exp(as.vector(rongelap_grid_pred))
# 线性预测的方差
rongelap_grid_var <- get_predVar(fit_rongelap_gaussian,
newdata = rongelap_grid_df, variances = list(predVar = TRUE), which = "predVar"
)
rongelap_grid_df$var_sk <- as.vector(rongelap_grid_var)
# 泊松分布
# 预测
rongelap_grid_pred <- predict(fit_rongelap_poisson,
newdata = rongelap_grid_df, type = "response"
)
rongelap_grid_df$pred_sp <- as.vector(rongelap_grid_pred)
# 线性预测的方差
rongelap_grid_var <- get_predVar(fit_rongelap_poisson,
newdata = rongelap_grid_df, variances = list(predVar = TRUE), which = "predVar"
)
rongelap_grid_df$var_sp <- as.vector(rongelap_grid_var)核辐射强度分布：对数高斯模型
图 11: 朗格拉普岛核辐射强度的分布
核辐射强度分布：泊松模型
图 12: 朗格拉普岛核辐射强度的分布
总结
- 空间数据探索
- 描述响应变量的分布
- 响应变量的空间分布
- 响应变量与横纵坐标的关系
- 建立广义线性模型以及检查模型残差分布
- 建立空间（广义）线性混合效应模型
- 模型拟合、评估和结果解释
- 准备未采样的位置数据
- 模型预测及展示预测结果
参考文献
Diggle, P. J., J. A. Tawn, 和 R. A. Moyeed. 1998. 《Model-based geostatistics》. Journal of the Royal Statistical Society: Series C (Applied Statistics) 47 (3): 299–350. https://doi.org/10.1111/1467-9876.00113.
Pebesma, Edzer. 2018. 《Simple Features for R: Standardized Support for Spatial Vector Data》. The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
Pebesma, Edzer, 和 Roger Bivand. 2023. Spatial Data Science: With applications in R. London: Chapman; Hall/CRC. https://doi.org/10.1201/9780429459016.
Rousset, François, 和 Jean-Baptiste Ferdy. 2014. 《Testing environmental and genetic effects in the presence of spatial autocorrelation》. Ecography 37 (8): 781–90. https://dx.doi.org/10.1111/ecog.00566.
谢谢
预测核辐射强度的分布 空间数据分析 黄湘云 统计之都 @ 北京会议中心 2023年7月12日