用 GRACE RL06 Level-2 月解和 GLDAS 水文模型数据反演地表质量异常和陆地水储量变化


以三峡水库(Three Gorges Reservoir, TGR)为例,用 2003 年 1 月份到 2016 年 8 月份的重力恢复及气候实验卫星(Gravity Recovery and Climate Experiment, GRACE) RL06 Level-2 月解以及全球陆地数据同化系统(Global Land Data Assimilation System, GLDAS)的水文模型数据(降水、蒸散和径流),结合退卷积法(Deconvolution)、水平衡方程(Water Balance Equation, WBE)与高斯过程回归法(Gaussian Process Regression, GPR)反演三峡水库及其周边区域的地表质量异常变化(Change of the Surface Mass Anomaly)和陆地水储量变化(Change of the Terrestrial Water Storage)。

数据准备

GRACE RL06 月解与 SLR RL06 月解可在地球科学信息系统数据中心(Information System and Data Center, ISDC) 或 NASA 的 EARTHDATA 下载。GLDAS 数据可从 GES DISC(Goddard Earth Sciences Data and Information Services Center) 下载,详情请参见 GLDAS 数据下载及使用。全球地形起伏模型 ETOPO5 可在美国国家海洋大气局(National Oceanic and Atmospheric Administration, NOAA) 下载或点击这里进行下载。全球陆地植被覆盖模型(Global Land Cover Map for 2009) 可在欧空局(European Space Agency, ESA)下载。

GRACE 数据处理

读取 SLR RL06 月解数据。注意:SLR RL06 月解需要截取到与 GRACE RL06 月解相对应的时间段,即从 2003.000 开始到 2016.5832 结束,对应 2003 年 1 月份到 2016 年 8月份,共 149 个 $C_{20}$ 系数。详情请参见地球引力场文件的读取与处理

读取 UT-CSR、GFZ、JPL 三家机构发布的 GRACE RL06 月解(见下表),然后用 SLR RL06 月解替换 GRACE RL06 月解中的 $C_{20}$,最后扣除球谐系数的平均值,得到去平均化的球谐系数,也即引力势异常的球谐系数。对三家去平均化的球谐系数求算数平均值,得到最终的球谐系数。

数据发布中心版本(Level-2)时间跨度球谐阶数参考引力场模型永久潮
UT-CSRRL062003 年 1 月 $\sim$ 2016 年 8 月96GGM05C
GFZRL062003 年 1 月 $\sim$ 2014 年 11 月60EIGEN-6C4不含
JPLRL062003 年 1 月 $\sim$ 2016 年 8 月60GGM05C

DDK5 去相关滤波处理引力势异常以光滑条带噪声。

将引力势异常转换到以等效水厚(EWT)表示的地表质量异常(SMA),详情请参见等效水高的计算和可视化

生成选定区域内的地表质量异常(SMA)。

点扩散函数对地表质量异常作退卷积处理,将扩散的地表质量异常恢复到原始的地表质量异常。当然,地表质量异常的恢复目前有两种比较成熟的方法:尺度因子法和 Mascon 方法。尺度因子方法对水文模型的模式有较大的依赖性,而 Mascon 方法的计算量偏大。退卷积法避免了地表质量异常对水文模型的依赖,同时计算量较小。

线性拟合(迭代拟合以剔除 RMS 大于 $3\sigma$ 的野值)选定区域内的地表质量异常,得到其变化率。

可视化选定区域的地表质量变化率。

GLDAS 数据处理

GLDAS 水文模型数据的时间范围为 2003 年 1 月份到 2016 年 8 月份,共 164 个月文件;空间范围为 $179.5^{\circ}W \sim 179.5^{\circ}E$,$59.5^{\circ}S \sim 89.5^{\circ}N$,分辨率为 $1^{\circ}\times 1^{\circ}$,网格的维度为 $150\times 360$。陆地水储量变化(TWSC)可由水平衡方程来估算,其表达式为
$$H_k = \frac{1}{\rho_w}\sum_{i=1}^{k} (P_i-ET_i-R_i)$$,其中 $H_k$ 为第 $k$ 个月相对于第 0 个月的陆地水储量变化(等效水厚);$k$ 的取值为 $1 \sim n$;$n$ 为总月数;$P_i$ 为第 $i$ 个月的降雨量与降雪量之和;$ET_i$ 为第 $i$ 个月的土壤水分蒸发量与植被蒸腾量之和;$R_i$ 为第 $i$ 个月的地表径流量与地下径流量之和。同计算地表质量异常的方法相似,对陆地水储量变化作去平均化处理。

修改陆地水储量变化的空间覆盖范围为 $0.5^{\circ}E \sim 359.5^{\circ}E$,$89.5^{\circ}N \sim 89.5^{\circ}S$,网格的维度为 $180\times 360$。

将陆地水储量变化插值到地表质量异常的网格,使陆地水储量变化的空间覆盖范围为 $0^{\circ}E \sim 359^{\circ}E$,$90^{\circ}N \sim 89^{\circ}S$,网格的维度为 $180\times 360$。

球谐展开陆地水储量变化到 89 阶,然后作高斯滤波光滑处理。需要说明的是 GLDAS 水文模型在海洋区域和南极大陆(约占全球表面积的 74%)不存在任何数据,简单地将 74% 的不存在的数据设为零会极大地缩减陆地数据的权重,因此这里采用了窗口球谐的方式展开陆地水储量变化,详情请参见球谐展开与谱分析。因为陆地水储量变化不含条带噪声,用 DDK5 去相关滤波对其光滑处理反而会添加条带噪声,因此这里用与 DDK5 光滑效果相当的滤波半径为 180km 的高斯滤波处理陆地水储量变化,以保持与 GRACE 数据处理的一致性。陆地窗口的构建如下:将全球地形起伏模型 ETOPO5 插值到 GLDAS 网格,然后将起伏值大于零的区域设置为陆地,其余的为海洋。由于陆地面积(排除南极大陆)约占全球面积的 26%,因此,对于展开到 89 阶的陆地窗口,需要 2106 个 Slepian 基函数来完成球谐恢复。

输出结果为

窗口球谐展开陆地水储量变化到 89 阶,然后用滤波半径为 180km 的高斯滤波对其作光滑处理。由于下面代码运行时间较长(大概需要 15 分钟),因此为了以后节约时间,我们将代码生成的“高斯滤波光滑后的陆地水储量变化”保存到文件  TWSCs_interp_grids_gaussian.npy 中,以方便下次使用。

加载 TWSCs_interp_grids_gaussian.npy 文件。

生成选定区域内的陆地水储量变化。

点扩散函数对陆地水储量变化作退卷积处理,将扩散的陆地水储量变化恢复到原始的陆地水储量变化。

线性拟合选定区域内的陆地水储量变化,得到其变化率。

可视化选定区域的陆地水储量变化率。

可视化选定区域 SMA-TWSC 的变化率,即地表质量异常变化率减去陆地水储量变化率。

根据 SMA-TWSC 变化率的分布来选择研究区域,例如包含三峡水库的一个多边形。

计算研究区域的面积。

输出结果为

计算研究区域内的地表质量异常。

线性拟合与高斯过程回归拟合研究区域内的地表质量异常时间序列,并可视化其变化趋势。

输出结果为

Model: GP regression
Objective: 531.8344574562797
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

GP_regression. valueconstraintspriors
rbf.variance 106.85016469358371 +ve
rbf.lengthscale 5.018094781335761 +ve
Gaussian_noise.variance 53.27703127219889 +ve

计算研究区域内的陆地水储量变化。

线性拟合与高斯过程回归拟合研究区域内的陆地水储量变化时间序列,并可视化其变化趋势。

输出结果为

Model: GP regression
Objective: 449.43313218184954
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

GP_regression. valueconstraintspriors
rbf.variance 23.998542855366136 +ve
rbf.lengthscale 1.7898924962292773 +ve
Gaussian_noise.variance 4.99652164566459 +ve

对陆地水储量的变化曲线作傅立叶频谱分析,得到其频谱分布。

研究区域内地表质量异常减去陆地水储量变化。

输出结果为

Model: GP regression
Objective: 529.1965792725174
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

GP_regression. valueconstraintspriors
rbf.variance 99.35388858495038 +ve
rbf.lengthscale 3.598655866056596 +ve
Gaussian_noise.variance46.466051027126916 +ve

可视化选定区域的地表质量变化率。

可视化选定区域陆地水储量变化率。

可视化选定区域 SMA-TWSC 的变化率。

三峡水库地形图和陆地植被覆盖图

三峡水库地形图
三峡水库植被覆盖

输出结果为

读取 Band

输出结果为

加载颜色表

输出结果为

读取 sheet

输出结果为

设置 color map,并作图。

输出图像为

由于受冰期回弹影响的区域主要集中在北美、斯堪的纳维亚、格陵兰以及南极等地区,因此这里没有对三峡水库及其周边区域进行冰川均衡调整(Glacial Isostatic Adjustment, GIA)改正;点扩散函数与纬度密切相关,在研究区域不大的情况下,可使用一个不变的点扩散函数;退卷积运算没有考虑球面的曲率效应,但这对局部区域质量异常的反演影响不大。


参考

致谢

根据邓梓锋同学的反馈,win 系统下安装 pyshtools 包(conda 环境下 pip install pyshtools) 显然有较多麻烦,原因可能是由于 pyshtools 是 fortran 程序 和 python 程序的结合体,若系统没有 fortran 语言编译器(gfortran、f95 等),可能会出现错误提示。针对 win 系统下安装 pyshtools 困难的问题,邓同学在经过多次尝试和反复确认之后,得出以下结论:

  • win 系统下安装 Fortran 编译器无法让 pyshtools 成功安装;
  • pyshtools 需要在 Linux 平台才能稳定运行;
  • 推荐使用 VMware Workstation Pro 15.0 的虚拟机进行操作,虚拟机安装的系统推荐使用 Ubuntu 18.04 及以上的版本。该系统版本入门比较简单,VMware Workstation Pro对该系统版本的优化也做得比较好,且能较好地运行pycharm。

关于 “用 GRACE RL06 Level-2 月解和 GLDAS 水文模型数据反演地表质量异常和陆地水储量变化” 的 4 个意见

  1. 博主你好,请问一下你在GRACE数据处理中为什么没有考虑一阶项系数的改正呢,是不是RL06模型已经将其用相应模型进行替换了呢?

    1. 一阶项系数(C10、C11、S11)的变化直接反映了地球质心的运动,尽管 GRACE 数据中默认将这些系数设为零,但实际上这些系数不为零,因此应当对其进行修正。博文中没有考虑一阶项系数的改正,是以下原因导致的:首先,一直寻求对应 RL06 数据产品的一阶项系数改正,但除了 Cheng 等人发布的 SLR-RL06 产品中有地球质心偏置外,目前 Swenson 等人貌似还没有发布直接对应 RL06 产品的一阶项系数改正。早在 2018 年 6 月份,给 Chambers 等人发邮件询问升级 RL06 产品中一阶项系数改正的问题,收到以下回复 “We are working on it, but we first have to evaluate the new coefficients and the new geocenter update. I doubt there will be a major change, so in the meantime just use the RL05 geocenter solutions.” 第二,目前,Swenson 等人可能已经发布了直接对应 RL06 产品的一阶项系数改正,但国内下不到任何数据。第三,用 Cheng 等人发布的 SLR-RL06 中的地球质心偏置反推一阶项会高估一阶项系数。第四,一阶项系数反映的空间尺度为半个球面,对于小尺度的质量变化,一阶项系数的改正与否影响不大。第五,球谐展开大地水准面扰动并进行去平均化,一阶项系数 C10、C11、S11 的周年振幅分别为 1.9mm、1.9
      mm、2.5mm,而 C20 的周年振幅为 14.1mm,这说明一阶项系数的变化相对于二阶项较小。第六,用 Swenson 等人的方法计算对应 RL06 数据产品的一阶项系数改正恐怕需要一些时间,但就目前来看,现在不可行。最后,若是您有这方面的数据资料或者资源链接等,可否分享一下?

  2. 博士大哥, 在DDK过滤步骤中
    filt_SHC_geopotential_Cs[i,0] = filt_SHC_geopotential_C.cnm
    filt_SHC_geopotential_Cs[i,1] = filt_SHC_geopotential_C.snm
    我觉得你这里的意思是obvact转化之后的nparray格式,取过滤后的值吧,这两行代码运行不了,可以讲讲吗

    1. 将 filt_SHC_geopotential_Clm,filt_SHC_geopotential_Slm = octave.filterSH(Wbd,SHC_geopotential_Cs[i,0],SHC_geopotential_Cs[i,1]) 换成 filt_SHC_geopotential_Clm,filt_SHC_geopotential_Slm = octave.filterSH(Wbd,SHC_geopotential_Cs[i,0],SHC_geopotential_Cs[i,1],nout=2)

发表评论

电子邮件地址不会被公开。 必填项已用*标注