GRACE重力卫星缺失值插补及深度学习重构数据（二），缺失值插补方法对比及代码实现。该方法也适用于任何缺失数据的栅格或NC文件

最新推荐文章于 2024-10-14 00:15:00 发布

笑不语

最新推荐文章于 2024-10-14 00:15:00 发布

阅读量2.1k

点赞数 6

分类专栏：地信 R语言新星计划文章标签：深度学习重构机器学习

本文链接：https://blog.csdn.net/weixin_43935907/article/details/131265332

版权

R语言新星计划同时被 2 个专栏收录

55 篇文章

订阅专栏

地信

2 篇文章

订阅专栏

一、引言

GRACE（Gravity Recovery and Climate Experiment）卫星是美国国家航空航天局（NASA）和德国航空航天中心（DLR）合作运营的一组卫星，在2002年至2017年期间对地球引力场进行了高精度测量。由于月球的引力也会对地球引力场产生影响，因此GRACE卫星还被用来探测月球引力场的细微变化，以研究月球内部的构造和演化。

GRACE月时变迁模型用于描述月球引力场的时变特征，是研究月球结构和演化的重要工具。然而，由于各种原因（如仪器故障、数据传输中断等），GRACE卫星获取的数据中可能会存在缺失值。这些缺失值会影响月时变迁模型的精度和可靠性，因此有必要采取合适的方法进行填补。

本文将介绍一种常用的缺失值填补方法——三次多项式插值法，并探讨其在GRACE月时变迁模型中的应用。这个方法在MODIS卫星或者任何缺失数据的栅格文件，nc文件都适用。

二、三次多项式插值法简介

插值法是一种在给定数据点的前提下，通过构造描述这些数据点的函数来近似未知函数的方法。当数据点存在缺失值时，插值法可以通过其他已知数据点来填补缺失值，从而使得数据的连续性和完整性得到保证。

三次多项式插值法是插值法中较为常用的一种方法。其基本思想是，对于给定的数据点，通过构造一个三次多项式函数来近似未知函数，并满足函数在每个数据点处的函数值等于相应数据点的函数值。通过求解三次多项式的系数，可以得到插值函数，并用于填补数据缺失值。

三次多项式插值法的优点在于可以较好地拟合数据的非线性特征，同时又不会过度拟合数据。此外，三次多项式插值法还具有较高的计算效率和较低的误差。

三、数据集

数据预处理

library(ncdf4)
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 <- na.spline(sampled_data$lwe_thickness)

# 过滤出缺失的数据，查看数据是否有插入
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.734405e+00
      2:   0.125 -89.625 2002-06-01            NA        -2.077785e+00
      3:   0.125 -89.375 2002-06-01            NA        -2.421165e+00
      4:   0.125 -89.125 2002-06-01            NA        -2.764544e+00
      5:   0.125 -88.875 2002-06-01            NA        -3.107923e+00
     ---                                                              
1036796: 359.875  88.875 2002-06-01            NA        -1.765162e+05
1036797: 359.875  89.125 2002-06-01            NA        -1.765162e+05
1036798: 359.875  89.375 2002-06-01            NA        -1.765162e+05
1036799: 359.875  89.625 2002-06-01            NA        -1.765162e+05
1036800: 359.875  89.875 2002-06-01            NA        -1.765162e+05

如果不想追加一列lwe_thickness_filled，可以直接赋值给lwe_thickness。我追加一列的目的是方便对比。