信号泄漏是用 GRACE GSM Level-2 月解反演地表质量变化的最大误差来源。信号泄漏是针对研究区域而言的,分为内部信号泄漏出研究区域(Leakage-out)和外部信号泄漏入研究区域(Leakage-in)。信号泄漏改正指将泄漏出的信号补回,泄漏入的信号扣除。若外部信号距离研究区域较远,那么只需选取较大的研究区域即可防止内部信号泄漏出区域边界。若选取的研究区域不够大,则内部的信号很容易扩散出区域边界,这时需要将泄漏出的信号补回。目前,主流的信号泄漏改正方法包括尺度因子法、正向模型法、谱域法、空域法,非主流的方法包括退卷积法和逆滤波法。下面以反演藏东南冰川融化为例,比较以上几种信号泄漏改正方法的异同。

信号泄漏改正前

下载、读取、去平均处理 CSR、GFZ、JPL 发布的 GRACE GSM RL06 Level-2 月解(96 阶),然后计算它们的算数平均值,得到综合解

输出结果为

结果显示综合解的时间覆盖范围为 2002-04 ~ 2019-10,月份总数为 211,月解总数为 178, 共缺失 33 个月解。之所以要求算数平均,是因为平均解的标准差比单家机构发布的月解的标准差小。设 CSR、GFZ、JPL 月解的标准差分别为 ​​$\sigma_{csr}$、​​$\sigma_{gfz}$、​​$\sigma_{jpl}$​​,那么综合解的标准差则为​ $\frac{\sqrt{\sigma^2_{csr}+\sigma^2_{gfz}+\sigma^2_{jpl}}}{3}$​​​。若三家机构的月解精度相当,均约等于 $\sigma$​,那么综合解的标准差为 $\frac{\sigma}{\sqrt{3}}$​,这比单家的标准差小了 ​$42\%$,因此,平均解的标准差比单家机构发布的月解的标准差小,更多详情请参见不同机构发布的 GRACE RL06 Level-2 月解的比较

下载、读取、去平均处理 SLR RL06 C20 月解

输出结果为

可视化 SLR RL06 C20 的变化

输出结果为

结果表明 $C_{20}$ 的变化具有周年项,长期表现为下降的趋势,其线性变化率为 $-1.0363 \times 10^{-11}/yr$。

由于 Lageos 等测地卫星(轨道高度约 5800km)相对于 GRACE 及 GRACE-FO(轨道高度约为 500km) 处于更高的位置,因此,其对地球的低阶引力位系数更为敏感,尤其是对地球动力学扁率的变化最敏感,因而一般用 SLR 测距数据解算的 C20 替换 GRACE 距离变化率数据解算的 C20。

输出结果为

由于 GRACE 及 GRACE-FO 卫星的轨道倾角为 $89^{\circ}$​,因此,在中低阶引力位系数的解算中,观测约束主要来自南北方向的星间距离变化率,因而造成了南北方向的条带效应。条带效应反映了高阶引力位系数之间存在相关性,消除条带噪声一种方法为去相关光滑滤波。这里采用 DDK 去相关滤波来压制南北条带噪声,更多详情请参见 DDK 去相关光滑滤波

输出结果为

结果显示 SLR C20 也经历了 DDK5 去相关滤波的处理,这样做的理由是后期改正信号泄漏误差时不必单独处理 C20。尽管 ICGEM 发布的 GRACE 月解是未替换 SLR C20 的 DDK 滤波解,然而先替换 C20 后滤波先滤波后替换 SLR C20 对结果几乎没有影响,详情请参考 GGTOOLS 使用文档

计算地表质量异常的变化。这里的异常是相对于参考引力场而言的,即 2002-04 ~ 2019-10 的平均引力场。因为这段时期内的平均引力场为一定值,不随时间变化,因此地表质量异常的变化与地表质量变化等价。

输出结果为

因为研究区域为藏东南的冰川融化区,因此,设置作图区域为藏东南及其周边地区。

输出结果为

结果显示研究区域占全球的面积比 $\alpha$ 为 ​$0.1\%$。之所以要计算面积比,是因为下文中用谱域法改正信号泄漏误差需要用到窗口球谐,而窗口球谐的 tapers(Slepian 基函数)数目是通过 Shannon 数 $N_s = (L+1)^2\alpha$ 所决定的,其中 $L$ 为球谐系数的最大阶数。以 Shannon 数为界,数量较大的 tapers 会微略地提高球谐恢复的准确度,但会极大地增加计算量,数量较小的 tapers 会显著地降低球谐恢复的准确度,因此,实际使用的 Shannon 数最好比理论计算的 Shannon 数稍大。

用迭代最小二乘法估算球谐系数(地表质量变化)的变化率,并将其展开到网格。注意:根据 Driscoll and Healy (1994) 样本定理,网格的分辨率是由球谐系数的最大阶数所决定的,两者满足关系 $N = 2(L+1)$,其中 $N$ 为沿纬度方向的节点数。96 阶的球谐系数可展开到维度为 $194\times388$ 的网格,维度为 $180\times360$ 的网格可通过插值得到。

输出结果为

计算研究区域内总地表质量变化的时间序列

输出结果为

结果显示藏东南地区的质量变化率为 $-6.7 ± 0.5$Gt。图中 “+” 号为时间序列,虚线为线性变化趋势,绿色实线为长期变化(线性变化加年际变化),红色实线为长期变化加周期变化。红色实线由高斯过程回归拟合得来,默认使用的核函数为 mat32+period。使用高斯过程回归的一个好处是不需要插值缺失的月解,尤其是在 2018 年前后 GRACE 与 GRACE-FO 交接时段有近一年的数据为缺失状态。

信号泄露改正

尺度因子法

基本原理:利用陆地水储量变化在高斯滤波前后的比例,恢复”地表质量变化”。”真地表质量变化”指假设高阶引力位系数之间不存在相关性,即没有条带噪声,无需进行去相关滤波处理。与之相对的是”地表质量变化”,即经过去条带噪声后的地表质量变化。

基本步骤:

  • 将 GLDAS 计算的陆地水储量展开到与 GRACE 同阶(例如 96 阶)的球谐系数
  • 对球谐系数进行高斯滤波处理,滤波半径需要与前面的 DDK 去相关滤波匹配
  • 计算研究区域内陆地水储量变化的时间序列
  • 计算尺度因子 $k = \mathbf y(t)/ \mathbf y_{gau}(t)$,其中 $ \mathbf y(t)$ 为高斯滤波前的时间序列,$ \mathbf y_{gau}(t)$ 为高斯滤波后的时间序列;可采用最小二乘法来拟合 $k$,即 $k = \frac{\mathbf y^T_{gau} \mathbf y}{\mathbf y^T_{gau} \mathbf y_{gau}}$

尺度因子法只有在满足 GRACE 反演的地表质量变化分布与 GLDAS 计算的陆地水储量变化分布相似的前提下才可以使用,否则尺度因子是不准确的,甚至是错误的。

尺度因子法示意图:上排图为研究区域(黑色实线为边界)内真质量分布(红色圈圈);下排图为滤波光滑后的结果

如上图所示,真质量分布经滤波后会发生扩散。若真质量分布远离研究区域的边界,那么扩散的质量分布一般不会超出研究区域(左图),此时的尺度因子等于 1;若真质量分布靠近边界,那么扩散后的质量分布会有一部分跨出边界(中图),此时的尺度因子大于 1;若边界外存在临近的质量,那么其扩散后会有一部分质量进入边界内,此时的尺度因子小于 1;以上仅为示例,只分析了同号质量分布情况,对于异号质量分布情况则更为复杂。总之,尺度因子与研究区域内外的质量分布密切相关,不同的质量分布会得到不同的尺度因子,具有明显的局部特性,并且依赖于水文模型。实际上,很难保证 GRACE 反演的地表质量变化分布与 GLDAS 计算的陆地水储量变化分布相似,因为地表质量变化不仅仅是陆地水储量,还包括地下水、地表水体(湖泊、水库)、冰川融化、泥沙淤积、冰期均衡调整、地震位错、以及其他地球物理过程引起的质量变化,而这些都是水文模型所不具备的,从而导致了地表质量变化分布与陆地水储量变化分布很难相似,除非研究区域内不存在这些因素。

输出结果为

结果显示尺度因子小于 1 且接近于 1,这说明滤波后有负的陆地水储量流出研究区域或有正的陆地水储量流入研究区域,具体是哪一种情况,需要作图进行判别。另外,该尺度因子到底能不能使用还需要看地表质量变化分布与陆地水储量变化分布是否相似,很遗憾两者并不相似,因此该尺度因子是不能使用的,也即地表质量变化率为 $-6.6$ Gt 是不准确的。

正向模型法

基本原理:”地表质量变化” $X$ 与”地表质量变化” $Y$ 之间可用函数 $Y = F(X)$ 来表达,其中 $F$ 代表高斯滤波。恢复”真地表质量变化”可理解为已知 $Y$ 和 $F$,求解 $X$。解算的方法有多种,在未知 $F$ 表达式的情况下,可采用迭代的方式解算,此即为正向模型法,又称为正向迭代法;在已知 $F$ 表达式的情况下,可采用反函数的方式解算,由于 $F$ 可表达为卷积的形式,因此该方法又称为退卷积法(见下文退卷积法)。

基本步骤:

  • 设地表质量变化(网格数据)的初始值为 $X_0$,然后将其展开到球谐系数
  • 计算 $F(X_0)$,设置 0 阶 和 1 阶球谐系数为 0(地球总质量不变,质心始终位于坐标系原点),然后将其展开到网格
  • 计算 $\Delta Y = Y – F(X_0)$,添加 $k\Delta Y$ 到 $X_0$,以更新初始值;$k$ 为增益因子,一般取值为 1.2;较大的值对应更快的收敛速度,但也降低了收敛的精度
  • 重复以上步骤,直到 $|\Delta Y|_{max}(\theta,\phi)$ 小于设定的限值,例如 10mm

输出结果为

谱域法

谱域法需要预先定义质量集中块(Mascon)的分布。对于冰川融化导致的地表质量变化,基于冰川目录选取面积大于 $10km^2$ 的冰川(glaciers.txt),然后根据这些冰川的位置生成质量块分布。这里选取面积大于 $10km^2$ 的冰川是因为面积小于$10km^2$ 的冰川固然也能生成 Mascon,但是这将会面临着由于 Mascon 数目较多而导致无法正确解算 Mascon 拟合值的问题。另一个原因是为了尽可能地保留冰川,生成尽可能多的质量块,让更多的冰川参与到地表质量拟合中。

基本原理:假设”真地表质量变化”是由一系列位置已知、质量未知的质量集中块(Mascon)引起的,根据质量块与球谐谱的线性关系(单位质量块对应谱 $f$,则值为 $m$ 的质量块对应谱 $mf$),已知”视地表质量变化”的谱可反推每一个质量块的值。

基本步骤:

  • 球谐展开单位质量块(例如边长为 $0.5^{\circ}$、等效水厚为 1mm、中心位于 $[\theta_i,\phi_i]$ 的 Mascon)并截断到特定阶数
  • 作高斯滤波处理,得到谱(球谐)系数 $f_i = [C_{0,0},…,C_{l,l}, S_{0,0},…,S_{l,l}]^T_i$
  • 若研究区域内存在 n 个质量块,分别记为 $m_1(\theta_1,\phi_1), …, m_n(\theta_n,\phi_n)$,谱系数为 $f_1, …, f_n$,则 $m_1f_1 + … + m_nf_n$ 为总的谱系数,矩阵形式为 $Am = b_{window}$,其中 $A = [f_1,…,f_l]$,$m = [m_1,…,m_n]^T$,$b_{window}$ 为窗口球谐处理后的 $b$,其中 $b = [C_{0,0},…,C_{l,l}, S_{0,0},…,S_{l,l}]^T$,为”视地表质量变化”的谱
  • 计算研究区域的 Slepian 基函数,然后对 $b$ 作窗口球谐处理,得到 $b_{wondow}$
  • 用  Tikhonov 正规化(病态的最小二乘估计)拟合参数 m

作窗口球谐处理是因”视地表质量变化”的谱是全球 Mascon 的综合结果,不仅仅是由研究区域内的 Mascon 引起的,还包含了研究区域以外的 Mascon。为了排除研究区域外的 Mascon 在”视地表质量变化”谱的分量,因此需要对谱作窗口处理。 用 Tikhonov 正规化拟合各个质量块的值,而不用普通的最小二乘法是因为 Tikhonov 正规化放宽了约束条件,能够避免 Mascon 的拟合值出现正负交错的现象,使得拟合值均为同号,尤其是在 Mascon 数目较多的情况下。

输出结果为

这里的标准差实际上不应为 0,只是因为窗口球谐的误差传播较为复杂,且暂时搁置了 Tikhonov 正规化中拟合参数的不确定度估计,目前待改进。

空域法

空域法与谱域法类似,也需要事先定义质量集中块的分布。

基本原理:假设”真地表质量变化”是由一系列位置已知、质量未知的质量集中块(Mascon)引起的,根据质量块与球谐展开到空间分布的线性关系(单位质量块对应谱的空间分布为 $f$,则值为 $m$ 的质量块对应 $mf$),已知”视地表质量变化”可反推每一个质量块的值。

基本步骤:

  • 球谐展开单位质量块(例如边长为 $0.5^{\circ}$、等效水厚为 1mm、中心位于 $[\theta_i,\phi_i]$ 的 Mascon)并截断到特定阶数
  • 作高斯滤波处理,得到谱(球谐)系数 $f_i = [C_{0,0},…,C_{l,l}, S_{0,0},…,S_{l,l}]^T_i$,然后展开到网格并保留研究区域内的网格 $g_i$
  • 若研究区域内存在 n 个质量块,分别记为 $m_1(\theta_1,\phi_1), …, m_n(\theta_n,\phi_n)$,球谐展开的空间分布为 $g_1, …, g_n$,则 $m_1g_1 + … + m_ng_n$ 为总的空间分布,矩阵形式为 $Am = b$,其中 $A = [g_1,…,g_l]$,$m = [m_1,…,m_n]^T$,$b$ 为研究区域内”视地表质量变化”的空间分布
  • 用  Tikhonov 正规化(病态的最小二乘估计)拟合参数 m

输出结果为

这里的标准差实际上不应为 0,只是因为暂时搁置了 Tikhonov 正规化中拟合参数的不确定度估计,目前待改进。

退卷积法

退卷积法可认为是正向模型法的另一支(见上文的正向模型法)。”真地表质量变化”经高斯滤波后会发生扩散(类似于光学中的点扩散函数或核函数),但扩散后的总量基本等于扩散前的总量。根据卷积定理,谱域上的积等于空域上的卷积,因此 $F(X) = X \otimes PSF$,其中 $PSF$ 为高斯滤波的逆傅立叶变换,也即点扩散函数。逆向思考,$X$ 等于滤波后的质量分布 $Y$ 与 $PSF$ 的退卷积。退卷积并不难实现,可采用 unsupervised Wiener-Hunt 算法计算退卷积,并且 scikit-image 程序包已经提供了该算法的实现。

逆滤波法

逆滤波法的思想来源于正向模型法和退卷积法。

基本原理:既然”视地表质量变化”是由”真地表质量变化”通过高斯滤波得来,而在谱域上,”视地表质量变化”的球谐系数等于高斯滤波系数乘以”真地表质量变化”的球谐系数,那么”真地表质量变化”的球谐系数只需要用”视地表质量变化”的球谐系数除以高斯滤波系数即可。

基本步骤:

  • “视地表质量变化”的球谐系数除以高斯滤波系数

输出结果为

从结果上看,逆滤波法与正向模型法符合地非常好。

结果的比较

到底哪一种方法更能反映实际情况,需要用实测的冰川质量平衡数据作验证(under construction),事实上,用谱域法计算的结果与实测数据很接近。用尺度因子法计算的结果偏小,并且在实际使用中会受到水文模型以及地表质量变化与陆地水储量变化分布是否相似的限制。空域法与谱域法是同一个事物的两面,两者都必须预先设定质量块的分布,且都使用了 Tikhonov 正规化。这两种方法的优点在于能直接反演是哪些质量块发生了变化。 与谱域法相比,空域法避免了窗口球谐的计算,节省了计算时间,但是计算结果偏大。正向模型法与逆滤波法的计算结果都是偏大。相对于正向模型法和退卷积法,逆滤波法及其简单,但效果并不输于两者。以上几种方法都不可避免地使用了高斯滤波,然而滤波半径的选取是至关重要的,这直接影响到信号泄漏改正。注意:退卷积法和逆滤波法均为非主流的方法,使用时请慎重。

研究区域内的冰川融化近似等于地表质量变化减去陆地水储量变化。除了谱域法和空域法外,其他的信号泄漏改正方法均可以先计算地表质量变化,再进行信号泄漏改正,然后扣除陆地水储量变化。由于谱域法和空域法默认质量块的变化全部由冰川融化导致,因此需要在 DDK 滤波前扣除陆地水储量变化。

19
Leave a Reply

avatar
7 Comment threads
12 Thread replies
0 Followers
 
Most reacted comment
Hottest comment thread
6 Comment authors
JIAOjhlDing XinJoanneKeepFi Recent comment authors
  Subscribe  
newest oldest most voted
Notify of
Jhualei
Guest
Jhualei

博士大哥你好!昨天无意中在网上看到你这个网站分享的东西很是感动,对我的帮助很大。我今天再浏览博文的时候发现排版混乱了,所有东西都是一行一行的,用同学的浏览器打开就可以正常浏览,我把网页的各种权限都打开还是不行,换火狐浏览器也不行。不知道你碰到过这个情况没有?

Jhualei
Guest
Jhualei

有点小激动,刚才还没把话说完,就是你的“信号泄漏误差修正方法的比较”这一篇文章我点击“继续阅读”是没有东西出来的,无法浏览,同学的试了一下也是这种情况,只显示“under construction”,这篇文章是你现在在更新维护吗?

Joanne
Guest
Joanne

博主您好,在看readme文件中,我使用
gsm_twsc_series = gsm_r_ddk5_sma_series_leakage – twsc_grid_series
gsm_twsc_series_rate = gsm_twsc_series.rate()

fig_name = ‘GRACE_TWSC_series.png’
ylabel = ‘GRACE-TWSC [Gt]’
gsm_twsc_series.plot(fig_name,ylabel,kernel=’rbf’)

出现了这种错误:
Traceback (most recent call last):
File “/home/han/PycharmProjects/GRACEdaima2ban-average (copy)/GRACEprocessing.py”, line 88, in
gsm_twsc_series = gsm_r_ddk5_sma_series_leakage – twsc_grid_series
File “/home/han/PycharmProjects/GRACEdaima2ban-average (copy)/class_Series.py”, line 342, in __sub__
if (self_existing_solution_flag != other_existing_solution_flag).any():
AttributeError: ‘bool’ object has no attribute ‘any’
在想了很久也在网上查了很久,实在不知道怎么解决,期待您的帮助,感谢博主。

KeepFi
Guest
KeepFi

博主您好,为什么下载数据服务器总是不响应呢,试了很多次不成功

Ding Xin
Guest
Ding Xin

博士,您好!我想问一下:为什么真地表质量变化经过高斯滤波得到视地表质量异常,有论文说为了减小高阶引力位系数误差,也需要使用高斯滤波。

jhl
Guest
jhl

李大哥春节好 我想问您一下求尺度因子k的那个小公式是怎么推导出来的 只知道利用最小二乘但是不知道怎么推导

JIAO
Guest
JIAO

描述的很简介漂亮. 谱域法对应于Jacob,2012? 空域法对应于Yi,2014? 不知理解的是否对, 或者是否有更好的文献参考可供推荐参考呢.thx 😉