如何在多图层Raster对象中逐像元去除时间序列异常值？

我有一个极高空间分辨率的NetCDF文件，参数如下：

dimensions : 2160, 4320, 9331200, 172  (nrow, ncol, ncell, nlayers)
resolution : 0.0833333, 0.0833333  (x, y)
extent     : -180, 179.9998, -90, 89.99993  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 

我的最终目标是计算全部172个图层（即时间维度）的均值。在此之前，我希望逐像元对时间序列去除异常值，例如使用公式：x ± 2.5 * std(x)，其中x为每个像元的时间序列数据。

我尝试将Raster对象转换为数据框，但由于分辨率过高，处理速度极慢。我也尝试过使用raster包中的clamp和reclassify函数，但这些函数使用的是针对整个Raster对象的固定阈值：

x <- clamp(x, lower = 0, upper = 2.5)

我需要的是基于每个像元时间序列（即上述动态上下限）去除异常值的函数。

是否可以直接在Raster格式下实现？无论是使用Raster包还是Terra包？

内容的提问来源于stack exchange，提问作者Francesco Giardina