以三峡水库(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)改正;点扩散函数与纬度密切相关,在研究区域不大的情况下,可使用一个不变的点扩散函数;退卷积运算没有考虑球面的曲率效应,但这对局部区域质量异常的反演影响不大。


5 1 vote
Article Rating
Subscribe
Notify of
guest
35 Comments
oldest
newest most voted
Inline Feedbacks
View all comments
廖洛
廖洛
April 16, 2019 17:07

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

han
han
November 15, 2019 20:35

博士大哥, 在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格式,取过滤后的值吧,这两行代码运行不了,可以讲讲吗

Chloe
Chloe
November 26, 2019 17:31

老师你好,
按照你写的步骤进行计算,在DDK滤波后得到EWT结果中第2015.02月的值明显比其他时期的值高/低很多,例如:
2015.01 范围:-2238~1105
2015.02 范围:-4235~6176
想请问一下,老师或者其他小伙伴有没有遇到这个问题。使用GLDAS计算的TWS没有这样的情况,不知道这种情况是计算误差还是别的原因导致的。

刘斌
刘斌
December 4, 2019 18:51

老师您好,GRACE球谐系数以及GLDAS计算的陆地水变化,您都做了去平均处理,请问为什么经过这一步处理?

Viva La Vida
Viva La Vida
December 10, 2019 18:59

老师您好,在
# Shannon Number
N = 2106 # (89+1)**2*0.26 = 2106

land_windows = SHWindow.from_mask(land_mask,lmax,N)
print(‘land windows’)
land_windows.info()
中,land_windows = SHWindow.from_mask(land_mask,lmax,N)出现了pyshtools.shtools.SHToolsError:Error allocating memory.
Error — SHReturnTapersMap
Problem allocating arrays DIJ and EVEC
请问是什么原因以及怎么处理,感谢老师

DCA
DCA
December 20, 2019 09:47

博主您好,怎么把研究区域换成自己想要的区域呢?

wuttingtao
wuttingtao
December 20, 2019 14:17

李老师你好,在RL06数据读取时,我用的是完整的三家机构的RL06 60阶的数据,数据的时长是一样,前边的代码已经针对性的修改过了,是不是就不需要下边贴的这部分代码了?其中[-16:,:,61:,61:]这个范围里的-16:,一直没看懂是什么意思,求解答。。。
SHC_geopotential_Cs[:-16,:,61:,61:] = SHC_geopotential_Cs[:-16,:,61:,61:]*3
SHC_geopotential_Cs[-16:,:,61:,61:] = SHC_geopotential_Cs[-16:,:,61:,61:]*2

SHC_geopotential_Cs_std[:-16,:,61:,61:]=SHC_geopotential_Cs[:-16,:,61:,61:]*3
SHC_geopotential_Cs_std[-16:,:,61:,61:]=SHC_geopotential_Cs[-16:,:,61:,61:]*2

wuttingtao
wuttingtao
December 20, 2019 16:53

李老师 你好,在文中多次出现了display(model),但未定义display
请问这里的display是从import librosa.dsiplay as display得到的吗?还是别的包里的函数。
感谢!

Rainbow
Rainbow
December 23, 2019 10:02

博主您好,我想了解一下,您反演的水储量空间分辨率是多少?展开到网格是否可以自己设定0.5°或者1°呢

taomy
taomy
January 15, 2020 19:59

博士,你好!在经过一段时间对您网站上的这篇文章学习后,有以下几个问题需要请教您,希望你能抽空回答。
(1)很多论文上说,GRACE数据经过DDK滤波后紧接着用半径为400-500km的高斯滤波去除高频噪声,为何您的文章上并没有采用高斯滤波处理?如果采用400-500km的高斯滤波处理,点扩散函数中高斯滤波半径应该如何选取?(是400-500km还是与DDK滤波对应的等效半径)
(2)“用与 DDK5 光滑效果相当的滤波半径为 180km 的高斯滤波处理陆地水储量变化”是不是为了去除噪声?滤波半径的选取仅仅是为了保持与GRACE的一致性吗?可否选取效果更好的滤波半径,比如上述的400-500km。
(3)对于研究区域,例如我输入的是某省的边界经纬度,是不规则的多边形,计算研究区域面积时出现变小的情况,这是计算面积出现错误还是把研究区域拟合成椭圆导致的呢?研究区域是否发生了变化?
(4)GLDAS中陆地水储量变化的傅立叶频谱分析的作用是什么?
(5)对于最终的计算结果SMA-TWSC,SMA-TWSC的时间序列(可视化时的‘’+‘’)与SMA-TWSC的变化趋势(可视化时的红线表示)两者选取谁作为最终结果。或者说,如果我有另外一份参考数据,我应该选取哪个作为我的最终计算结果与参考数据对比呢?(对了,博士知道官方发布的全国各省地下水变化的数据相应网址吗?)

Stone
Stone
March 28, 2020 19:14

老师您好,请问psf = np.load(‘intermediate/psf_tgr.npy’),这里的文件在哪里可以下载呢?或者怎么制作这个文件呢?谢谢

Stone
Stone
April 2, 2020 12:51

李老师您好,我在“用 DDK5 去相关滤波处理引力势异常以光滑条带噪声”代码
for i in range(SHC_geopotential_Cs.shape[0]):
filt_SHC_geopotential_C = octave.filterSH(Wbd,SHC_geopotential_Cs[i,0],SHC_geopotential_Cs[i,1],nout=2)
filt_SHC_geopotential_Cs[i,0] = filt_SHC_geopotential_C.cnm
filt_SHC_geopotential_Cs[i,1] = filt_SHC_geopotential_C.snm

中,请问.cnm和.snm函数有什么作用,是取filt_SHC_geopotential_C的哪部分赋给filt_SHC_geopotential_Cs?
我运行时出现以下错误:
filt_SHC_geopotential_Cs[i,0] = filt_SHC_geopotential_C.cnm
AttributeError: 'list' object has no attribute 'cnm'

我已经得到了filt_SHC_geopotential_C,能不能换个方法赋给filt_SHC_geopotential_Cs呢?
感谢老师

wkk
wkk
July 4, 2020 18:31

老师,您好,我在进行滤波处理时,显示如下报错信息,希望老师抽空看一下,感谢老师

777.PNG