一、引言
GRACE（Gravity Recovery and Climate Experiment）卫星是美国国家航空航天局（NASA）和德国航空航天中心（DLR）合作运营的一组卫星，在2002年至2017年期间对地球引力场进行了高精度测量。由于月球的引力也会对地球引力场产生影响，因此GRACE卫星还被用来探测月球引力场的细微变化，以研究月球内部的构造和演化。
GRACE月时变迁模型用于描述月球引力场的时变特征，是研究月球结构和演化的重要工具。然而，由于各种原因（如仪器故障、数据传输中断等），GRACE卫星获取的数据中可能会存在缺失值。这些缺失值会影响月时变迁模型的精度和可靠性，因此有必要采取合适的方法进行填补。
本文将介绍一种常用的缺失值填补方法——三次多项式插值法，并探讨其在GRACE月时变迁模型中的应用。这个方法在MODIS卫星或者任何缺失数据的栅格文件，nc文件都适用。
二、三次多项式插值法简介
插值法是一种在给定数据点的前提下，通过构造描述这些数据点的函数来近似未知函数的方法。当数据点存在缺失值时，插值法可以通过其他已知数据点来填补缺失值，从而使得数据的连续性和完整性得到保证。
三次多项式插值法是插值法中较为常用的一种方法。其基本思想是，对于给定的数据点，通过构造一个三次多项式函数来近似未知函数，并满足函数在每个数据点处的函数值等于相应数据点的函数值。通过求解三次多项式的系数，可以得到插值函数，并用于填补数据缺失值。
三次多项式插值法的优点在于可以较好地拟合数据的非线性特征，同时又不会过度拟合数据。此外，三次多项式插值法还具有较高的计算效率和较低的误差。
三、数据集
- 数据预处理
library(ncdf4) nstall.packages("rugarch") install.packages("zoo") library(rugarch) library(zoo) ncdata <- nc_open("data.nc") lon <- ncvar_get(ncdata,'lon') # 经度 lat <- ncvar_get(ncdata,'lat') # 纬度 time <- ncvar_get(ncdata,'time') # 时间 lwe_thickness <- ncvar_get(ncdata,'lwe_thickness') # 覆盖厚度 # 这份测绘数据起始时间是2002.01.01，time的值表示过去的天数 time <- as.Date("2002-01-01") + time library(data.table) # 对数据进行重组和转置，以扁平化数据结构 lwe_thickness_flat <- as.data.table(lwe_thickness) # lon lwe_thickness_flat$V1 <- lon[lwe_thickness_flat$V1] # lat lwe_thickness_flat$V2 <- lat[lwe_thickness_flat$V2] # time lwe_thickness_flat$V3 <- time[lwe_thickness_flat$V3] # 将数据框lwe_thickness_flat转换为数据表对象 lwe_thickness_flat<-as.data.table(lwe_thickness_flat) # 设置列名 setnames(lwe_thickness_flat, c("V1","V2","V3","value"), c("lon","lat","time","lwe_thickness"))
- 获取缺失子集作为演示
library(dplyr) # 对数据进行日期筛选 sampled_data <- lwe_thickness_flat %>% filter(time < as.Date("2003-12-16"))
- 结果展示
lon lat time lwe_thickness 1: 0.125 -89.875 2002-04-18 -2.468668 2: 0.125 -89.875 2002-05-10 -1.569404 3: 0.125 -89.875 2002-08-16 -3.044039 4: 0.125 -89.875 2002-09-16 -3.201635 5: 0.125 -89.875 2002-10-16 -1.785586 --- 17625596: 359.875 89.875 2003-07-16 3.973456 17625597: 359.875 89.875 2003-08-16 3.146194 17625598: 359.875 89.875 2003-09-16 3.110446 17625599: 359.875 89.875 2003-10-16 4.867654 17625600: 359.875 89.875 2003-11-16 2.196712
四、填补缺失值
- 找出缺失的年月
time_sub <- unique(sampled_data$time) # 提取年月信息 year_month <- format(as.Date(time_sub), "%Y-%m") # 构造包含所有年月的数据框 start_date <- min(year_month) end_date <- max(year_month) all_dates <- seq(as.Date(paste0(start_date, "-01")), as.Date(paste0(end_date, "-01")), by = "month") all_year_month <- format(all_dates, "%Y-%m") df_all_dates <- data.frame(year_month = all_year_month, date = all_dates, stringsAsFactors = FALSE) # 找出缺失的日期 missing_dates <- df_all_dates[!df_all_dates$year_month %in% year_month, ] # 打印输出缺失的日期 print(missing_dates$year_month)
- 填充经度、纬度和时间
# 选取一个栅格的数据 grid_data <- lwe_thickness_flat %>% filter(time == as.Date("2002-04-18")) grid_data$lwe_thickness <- NA # 插入缺失年月的数据 for(i in 1:nrow(missing_dates)){ grid_data$time <- missing_dates$date[i] sampled_data <- rbind(sampled_data, grid_data) }
- 按照时间排序
sampled_data <- sampled_data[order(sampled_data$time)] # 检查排序结果 time_order <- unique(sampled_data$time) time_order
结果展示：
[1] "2002-04-18" "2002-05-10" "2002-06-01" "2002-07-01" "2002-08-16" [6] "2002-09-16" "2002-10-16" "2002-11-16" "2002-12-16" "2003-01-16" [11] "2003-02-15" "2003-03-16" "2003-04-16" "2003-05-11" "2003-06-01" [16] "2003-07-16" "2003-08-16" "2003-09-16" "2003-10-16" "2003-11-16"
缺失的三个月份补充进去了
- 三次多项式插值
# 导入插值包 install.packages("rugarch") install.packages("zoo") library(rugarch) library(zoo) library(imputeTS) # 使用三次多项式插值填补缺失值 sampled_data$lwe_thickness_filled <- approx(sampled_data$lwe_thickness, n = length(sampled_data$lwe_thickness))$y # 过滤出缺失的数据，查看数据是否有插入 grid_data_filled <- sampled_data %>% filter(time == as.Date("2002-06-01")) # 写入文件并忽略行名 fwrite(sampled_data, file = "D:/log/grace.csv", row.names = FALSE)
结果展示：
lon lat time lwe_thickness lwe_thickness_filled 1: 0.125 -89.875 2002-06-01 NA -1.391025 2: 0.125 -89.625 2002-06-01 NA -1.391026 3: 0.125 -89.375 2002-06-01 NA -1.391027 4: 0.125 -89.125 2002-06-01 NA -1.391028 5: 0.125 -88.875 2002-06-01 NA -1.391028 --- 1036796: 359.875 88.875 2002-06-01 NA -2.217528 1036797: 359.875 89.125 2002-06-01 NA -2.217529 1036798: 359.875 89.375 2002-06-01 NA -2.217530 1036799: 359.875 89.625 2002-06-01 NA -2.217531 1036800: 359.875 89.875 2002-06-01 NA -2.217531
如果不想追加一列lwe_thickness_filled，可以直接赋值给lwe_thickness。我追加一列的目的是方便对比。
五、总结
三次多项式插值使用的是拉格朗日插值法或牛顿插值法，这两种方法都是将给定数据点用多项式函数连接起来，从而进行插值。
对于拉格朗日插值法，要求使用的多项式为n次多项式，其中n为插值点的数量。其基本思路是：设待求函数为f(x)，通过构造n次多项式L(x)来逼近f(x)，满足在n个给定点xi处，L(x)与f(x)相等，且该多项式在所有点x上具有连续性和可导性。因此，通过L(x)来近似f(x)的值来实现插值。
对于牛顿插值法，同样也要求使用的多项式为n次多项式。它是通过将多项式表示成一系列插值点的差商来近似函数f(x)。与拉格朗日插值法不同，牛顿插值法在构造多项式时不需要解出系数，而是通过求导数来递归地表达所有的差商。
虽然两种方法的构造方式不同，但它们都可以被用来实现三次多项式插值。同时需要注意，当插值点数量增加时，这两种方法的计算量也会随之增加。
参考文献：
Using Satellite‑Based Terrestrial Water Storage Data: A Review
