在 R-Studio 中使用 gstat 进行时空插值 - 拟合正确的变异函数

如何解决在 R-Studio 中使用 gstat 进行时空插值 - 拟合正确的变异函数

我对空间评估非常陌生，来自心理学。

我正在使用软件 R 和包“gstat”和“spacetime”。

我想做一个时空插值。为此，我遵循 Gräler 等人的论文。 (https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf)

很遗憾，我找不到/适合正确的变异函数模型。我可以创建经验变异函数，这对我来说也是结论性的，但我没有进一步了解。我不明白如何定义各个参数，例如“sill”或“nugget”或它们代表什么。

这是我以前的方法：

我的 ST 数据框：

> 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,sill = 124)

我明白，我必须为空间和时间创建一个变异函数。但我不知道如何正确定义变异函数的参数（psill、模型、范围、金块）。如果我使用与论文中相同的参数，我会得到以下图：

如您所见，在此模型中，我只有两个滞后，变异函数看起来很奇怪。所以，我认为这是因为我的模型参数错误。我还测试了另一种方法，如论文中所述。

sumMetricModel <- vgmST("sumMetric",space = vgm(20,150,1),time = vgm(10,2,0.5),joint = vgm(80,1500,2.5),stAni = 120)

同样，我不知道如何设置参数，所以我采用了与论文中相同的方法并绘制了所有三个参数。

我很确定这取决于参数，但我真的不知道如何找到正确的变异函数模型来拟合。

第一个回答后更新：

正如建议的那样，我为新的经验变异函数使用了更少的时间间隔：

eVgmPm10<-variogramST(MMW10~1,tlags = 0:2)

我得到了这个线框：

如您所见，我从“滞后 0”开始有 3 个时间步长。

我尝试按照建议填充模型参数。

separableModel<-vgmST("separable",space=vgm(psill = 15,model = "Exp",range =250,nugget = 5),time = vgm(psill = 15,range =1,nugget = 0),sill = 15)

这是两者的线框图：

我真的很沮丧，也很抱歉我的知识不足，但我真的很努力去理解它，我很感激每一个帮助

解决方法

sill、range 和 nugget 的初始参数可以从经验时空变异函数中读取。通过以下方式生成经验变异函数的 3D 线框图：

plot(eVgmPm10,wireframe=T,scales=list(arrows=F))

金块是关于表面与 z 轴的近似交点（这相当于模型无法解释的小尺度变化）。门槛大约是您的经验曲面接近的最大 z 值，范围大约是曲面变平的距离。您将拥有空间（x 轴）和时间（y 轴）范围。

不同的模型对这些值的解释略有不同；详细信息在您在问题中引用的 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)

参考您问题的标题“正确的变异函数”，很难确定它是什么。您只能找到一个模型，该模型将是您提供给算法的数据中依赖结构的“良好”近似。例如，这种依赖性可能会在冬季和夏季月份之间发生变化，并且平均值中可能存在季节性成分，这些成分也需要反映在模型中。