预测核辐射强度的分布

目录

1. 核污染研究的背景和意义

2. 数据探索和分析的过程

   • 描述数据的空间分布
   • 广义线性模型及模型诊断
   • 空间广义线性混合效应模型
   • 模型拟合及准备预测位置
   • 模型预测及展示预测分布

氢弹核试验

朗格拉普岛

图 1: 朗格拉普岛地理位置特殊，是一个重要的军事基地，建有环岛公路和机场。卫星拍摄的影像数据（2020 年）来自 https://rmi-data.sprep.org/

采集的核辐射数据（部分）

采样数据的空间分布

岛上各采样点的核辐射强度

图 4: 岛上各采样点的核辐射强度

广义线性模型及拟合

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

空间广义线性混合效应模型

广义线性模型并没有考虑距离相关性，它认为各个观测点的数据是相互独立的。因此，在广义线性模型的基础上添加位置相关的随机效应，用以刻画未能直接观测到的潜在影响。根据 {}^{137}\mathrm{Cs} 放出伽马射线，在 n=157 个采样点，分别以时间间隔 t_i 测量辐射量 y(x_i)，建立泊松型空间广义线性混合效应模型(Diggle, Tawn, 和 Moyeed 1998)。

\begin{aligned} \log\{\lambda(x_i)\} & = \beta + S(x_{i})\\ y(x_{i}) &\sim \mathrm{Poisson}\big(t_i\lambda(x_i)\big) \end{aligned} \tag{1}

其中，\beta 表示截距，相当于平均水平，\lambda(x_i) 表示位置 x_i 处的辐射强度，S(x_{i}) 表示位置 x_i 处的空间效应，S(x),x \in \mathcal{D} \subset{\mathbb{R}^2} 是二维平稳空间高斯过程 \mathcal{S} 的具体实现。 \mathcal{D} 表示研究区域，可以理解为朗格拉普岛，它是二维实平面 \mathbb{R}^2 的子集。

随机过程的自协方差函数

随机过程 S(x) 的自协方差函数常用的有指数型、幂二次指数型（高斯型）和梅隆型，形式如下：

\begin{aligned} \mathsf{Cov}\{S(x_i), S(x_j)\} &= \sigma^2 \exp\big( -\frac{\|x_i -x_j\|_{2}}{\phi} \big) \\ \mathsf{Cov}\{ S(x_i), S(x_j) \} &= \sigma^2 \exp\big( -\frac{\|x_i -x_j\|_{2}^{2}}{2\phi^2} \big) \\ \mathsf{Cov}\{ S(x_i), S(x_j) \} &= \sigma^2 \frac{2^{1 - \nu}}{\Gamma(\nu)} \left(\sqrt{2\nu}\frac{\|x_i -x_j\|_{2}}{\phi}\right)^{\nu} K_{\nu}\left(\sqrt{2\nu}\frac{\|x_i -x_j\|_{2}}{\phi}\right) \\ K_{\nu}(x) &= \int_{0}^{\infty}\exp(-x \cosh t) \cosh (\nu t) \mathrm{dt} \end{aligned}

待估参数：代表方差的 \sigma^2 和代表范围的 \phi 。当 \nu = 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) 对应模型参数 \beta = 1.829 ，随机效应中的方差参数 lambda 对应模型参数 \sigma^2 = 0.3069 ，随机效应中相关参数 1.rho 的倒数对应模型范围参数 \phi = 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: 朗格拉普岛核辐射强度的分布

总结

1. 空间数据探索
   • 描述响应变量的分布
   • 响应变量的空间分布
   • 响应变量与横纵坐标的关系
2. 建立广义线性模型以及检查模型残差分布
3. 建立空间（广义）线性混合效应模型
4. 模型拟合、评估和结果解释
5. 准备未采样的位置数据
6. 模型预测及展示预测结果