在 R-Studio 中使用 gstat 进行时空插值 - 拟合正确的变异函数
Spatio Temporal Interpolation with gstat in R-Studio - Fitting the correct variogram
我对空间评估很陌生，来自心理学。
我正在使用软件 R 和软件包“gstat”和“spacetime”。
我想做一个时空插值。为此，我遵循 Gräler 等人的论文。 (https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf)
不幸的是，我不能 find/fit 正确的变差函数模型。我可以创建经验变异函数，这对我来说也是决定性的，但我没有进一步了解。我不明白如何定义单独的参数，例如“基台”或“金块”或它们代表什么。
这是我以前的方法：
我的 ST-Dataframe：
> df.stf
An object of class "STFDF"
Slot "data":
# A tibble: 2,406 x 6
KRS month DeprIndex MMW10 MMW25 NewInfect
<dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 1001 1 -4.08 NA NA NA
2 1002 1 -2.28 13.3 NA NA
3 1003 1 -3.29 17.8 11.7 NA
4 1004 1 -4.31 NA NA NA
5 1051 1 -2.51 17.7 NA NA
6 1053 1 -1.07 13.0 9.93 NA
7 1054 1 -0.863 NA NA NA
8 1055 1 -1.63 NA NA NA
9 1056 1 0.887 NA NA NA
10 1057 1 -1.21 14.7 8.84 NA
# ... with 2,396 more rows
Slot "sp":
class : SpatialPolygonsDataFrame
features : 401
extent : 3280359, 3921536, 5237511, 6103443 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs
variables : 1
names : AGS
min values : 1001
max values : 16077
Slot "time":
timeIndex
0001-01-30 1
0002-01-30 2
0003-01-30 3
0004-01-30 4
0005-01-30 5
0006-01-30 6
Slot "endTime":
[1] "0002-01-30 CET" "0003-01-30 CET" "0004-01-30 CET" "0005-01-30 CET" "0006-01-30 CET" "0007-01-30 CET"
所以，我有 401 个德国县和六个月内的月平均 PM10 值。
我的经验变异函数：
eVgmPm10<-variogramST(MMW10~1,df.stf,tlags = 0:5)
eVgmPm10$dist<-eVgmPm10$dist/1000
eVgmPm10$avgDist <- eVgmPm10$avgDist/1000
eVgmPm10$spacelag <- eVgmPm10$spacelag/1000
plot(eVgmPm10, map=F)
到目前为止一切正常（我猜）
现在我想像论文中解释的那样制作一个可分离的模型：
separableModel <- vgmST("separable",
space = vgm(0.9, "Exp", 200, 0.1),
time = vgm(0.9, "Sph", 3.5, 0.1),
sill = 124)
我明白，我必须为 space 和时间创建一个变差函数。但我不知道如何为变差函数正确定义参数（psill、模型、范围、块）。如果我使用与论文中相同的参数，我会得到以下图：
如您所见，在这个模型中，我只有两个滞后，变异函数看起来很奇怪。所以，我认为这是因为我的模型参数错误。我还测试了论文中的另一种方法。
sumMetricModel <- vgmST("sumMetric",
space = vgm(20, "Sph", 150, 1),
time = vgm(10, "Exp", 2, 0.5),
joint = vgm(80, "Sph", 1500, 2.5),
stAni = 120)
同样，我不知道如何设置参数，所以我采用了与论文中相同的方法并绘制了所有三个参数。
我很确定它取决于参数，但我真的不知道如何找到适合的正确变差函数模型。
第一个答案后更新：
按照建议，我为新的经验变差函数使用了更少的时间延迟：
eVgmPm10<-variogramST(MMW10~1,df.stf,tlags = 0:2)
我得到了这个线框图：
如您所见，我仍然有 3 个时间步从“滞后 0”开始。
然后我尝试按照建议填充模型参数。
separableModel<-vgmST("separable",
space=vgm(psill = 15,model = "Exp", range =250 ,nugget = 5),
time = vgm(psill = 15,model = "Exp", range =1 ,nugget = 0),
sill = 15)
这是两者的线框图：
我真的很沮丧，我很抱歉我的知识匮乏，但我真的很努力去理解它，我很感激每一个帮助
sill、range和nugget的初始参数可以从经验时空变异函数中读取。通过以下方式生成经验变异函数的 3D 线框图：
plot(eVgmPm10, wireframe=T, scales=list(arrows=F))
金块大约是表面与 z-axis 的近似交点（这相当于模型无法解释的小规模可变性）。基台大约是经验表面接近的最大值 z-value，范围大约是表面变平的距离。您将拥有空间 (x-axis) 和时间 (y-axis) 范围。
不同的模型对这些值的解释略有不同；详细信息在您在问题中引用的 gstat 小插图的第 2.1 节中给出。在可分离的情况下，例如，空间和时间变异函数分量被归一化为 1 的联合基台（基台 + 金块），并且模型被模型定义中的基台参数“拉伸”。在 sum-metric 模型中，此规范化不适用。对于度量分量，必须拟合时空各向异性，以便在单个模型中关联空间和时间距离（检查 gstat::estiStAni()）。我建议在 R 中使用 运行 demo("stkrige") 并试验这些参数以更好地了解它们的影响。请注意，参数也会相互影响，这有时会使解释和数值优化变得麻烦。
根据您在问题中提供的内容，5 的金块、15 的基台和 250 的空间范围可能是开始的值。此外，它看起来好像您的经验表面沿着时间域以波浪形式移动。这可能是数据集时间覆盖的人工制品。 6 个时间步长（在您的情况下为几个月）很少用于估计任何时间依赖性，尤其是超过 1 个（或可能是 2 个）步长。您的绘图表明您正在查看 5 个时间步长，其中“lag5”将仅基于第一个月和最后一个月的值。在估计经验变异函数时，请考虑使用较少的时间滞后：
eVgmPm10<-variogramST(MMW10~1, df.stf, tlags = 0:2)
关于您问题的标题“正确的变差函数”，很难确定那是什么。您只能找到一个模型，该模型将是您提供给算法的数据中依赖结构的“良好”近似。例如，依赖性可能会在冬季和夏季月份之间发生变化，并且平均值中可能存在季节性成分，这也需要反映在模型中。
我对空间评估很陌生，来自心理学。
我正在使用软件 R 和软件包“gstat”和“spacetime”。
我想做一个时空插值。为此，我遵循 Gräler 等人的论文。 (https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf)
不幸的是，我不能 find/fit 正确的变差函数模型。我可以创建经验变异函数，这对我来说也是决定性的，但我没有进一步了解。我不明白如何定义单独的参数，例如“基台”或“金块”或它们代表什么。
这是我以前的方法：
我的 ST-Dataframe：
> df.stf
An object of class "STFDF"
Slot "data":
# A tibble: 2,406 x 6
KRS month DeprIndex MMW10 MMW25 NewInfect
<dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 1001 1 -4.08 NA NA NA
2 1002 1 -2.28 13.3 NA NA
3 1003 1 -3.29 17.8 11.7 NA
4 1004 1 -4.31 NA NA NA
5 1051 1 -2.51 17.7 NA NA
6 1053 1 -1.07 13.0 9.93 NA
7 1054 1 -0.863 NA NA NA
8 1055 1 -1.63 NA NA NA
9 1056 1 0.887 NA NA NA
10 1057 1 -1.21 14.7 8.84 NA
# ... with 2,396 more rows
Slot "sp":
class : SpatialPolygonsDataFrame
features : 401
extent : 3280359, 3921536, 5237511, 6103443 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs
variables : 1
names : AGS
min values : 1001
max values : 16077
Slot "time":
timeIndex
0001-01-30 1
0002-01-30 2
0003-01-30 3
0004-01-30 4
0005-01-30 5
0006-01-30 6
Slot "endTime":
[1] "0002-01-30 CET" "0003-01-30 CET" "0004-01-30 CET" "0005-01-30 CET" "0006-01-30 CET" "0007-01-30 CET"
所以，我有 401 个德国县和六个月内的月平均 PM10 值。
我的经验变异函数：
eVgmPm10<-variogramST(MMW10~1,df.stf,tlags = 0:5)
eVgmPm10$dist<-eVgmPm10$dist/1000
eVgmPm10$avgDist <- eVgmPm10$avgDist/1000
eVgmPm10$spacelag <- eVgmPm10$spacelag/1000
plot(eVgmPm10, map=F)
到目前为止一切正常（我猜）
现在我想像论文中解释的那样制作一个可分离的模型：
separableModel <- vgmST("separable",
space = vgm(0.9, "Exp", 200, 0.1),
time = vgm(0.9, "Sph", 3.5, 0.1),
sill = 124)
我明白，我必须为 space 和时间创建一个变差函数。但我不知道如何为变差函数正确定义参数（psill、模型、范围、块）。如果我使用与论文中相同的参数，我会得到以下图：
如您所见，在这个模型中，我只有两个滞后，变异函数看起来很奇怪。所以，我认为这是因为我的模型参数错误。我还测试了论文中的另一种方法。
sumMetricModel <- vgmST("sumMetric",
space = vgm(20, "Sph", 150, 1),
time = vgm(10, "Exp", 2, 0.5),
joint = vgm(80, "Sph", 1500, 2.5),
stAni = 120)
同样，我不知道如何设置参数，所以我采用了与论文中相同的方法并绘制了所有三个参数。
我很确定它取决于参数，但我真的不知道如何找到适合的正确变差函数模型。
第一个答案后更新：
按照建议，我为新的经验变差函数使用了更少的时间延迟：
eVgmPm10<-variogramST(MMW10~1,df.stf,tlags = 0:2)
我得到了这个线框图：
如您所见，我仍然有 3 个时间步从“滞后 0”开始。
然后我尝试按照建议填充模型参数。
separableModel<-vgmST("separable",
space=vgm(psill = 15,model = "Exp", range =250 ,nugget = 5),
time = vgm(psill = 15,model = "Exp", range =1 ,nugget = 0),
sill = 15)
这是两者的线框图：
我真的很沮丧，我很抱歉我的知识匮乏，但我真的很努力去理解它，我很感激每一个帮助
sill、range和nugget的初始参数可以从经验时空变异函数中读取。通过以下方式生成经验变异函数的 3D 线框图：
plot(eVgmPm10, wireframe=T, scales=list(arrows=F))
金块大约是表面与 z-axis 的近似交点（这相当于模型无法解释的小规模可变性）。基台大约是经验表面接近的最大值 z-value，范围大约是表面变平的距离。您将拥有空间 (x-axis) 和时间 (y-axis) 范围。
不同的模型对这些值的解释略有不同；详细信息在您在问题中引用的 gstat 小插图的第 2.1 节中给出。在可分离的情况下，例如，空间和时间变异函数分量被归一化为 1 的联合基台（基台 + 金块），并且模型被模型定义中的基台参数“拉伸”。在 sum-metric 模型中，此规范化不适用。对于度量分量，必须拟合时空各向异性，以便在单个模型中关联空间和时间距离（检查 gstat::estiStAni()）。我建议在 R 中使用 运行 demo("stkrige") 并试验这些参数以更好地了解它们的影响。请注意，参数也会相互影响，这有时会使解释和数值优化变得麻烦。
根据您在问题中提供的内容，5 的金块、15 的基台和 250 的空间范围可能是开始的值。此外，它看起来好像您的经验表面沿着时间域以波浪形式移动。这可能是数据集时间覆盖的人工制品。 6 个时间步长（在您的情况下为几个月）很少用于估计任何时间依赖性，尤其是超过 1 个（或可能是 2 个）步长。您的绘图表明您正在查看 5 个时间步长，其中“lag5”将仅基于第一个月和最后一个月的值。在估计经验变异函数时，请考虑使用较少的时间滞后：
eVgmPm10<-variogramST(MMW10~1, df.stf, tlags = 0:2)
关于您问题的标题“正确的变差函数”，很难确定那是什么。您只能找到一个模型，该模型将是您提供给算法的数据中依赖结构的“良好”近似。例如，依赖性可能会在冬季和夏季月份之间发生变化，并且平均值中可能存在季节性成分，这也需要反映在模型中。