利用重力卫星监测全球陆地水储量变化时不同滤波算法的影响

崔立鲁，张 诚，何明睿

(成都大学 建筑与土木工程学院，四川 成都 610106)

0 引 言

随着全球经济快速发展和人口的快速增加，导致全球气候变暖，极端天气频发，水资源在时空分布上呈现出严重的不均匀现象，极大地影响了人类生存和发展[1].随着2002年4月重力恢复与气候实验(Gravity Recovery and Climate Explorer,GRACE)卫星正式提供服务，一种全新的对地观测手段出现在各国科学家面前.相对于传统监测技术手段，该手段可以提供低成本、持续性、大面积的观测数据，且这些观测数据能够反映出陆地水储量变化(Terrestrial Water Storage Change,TWSC)总的情况[2].因此，科学家们利用GRACE对卫星观测数据进行了大量的水文应用研究[3-5].

由于卫星轨道误差、卫星载荷误差、仪器测量误差以及重力场模型本身缺陷等，由GRACE时变重力场模型所得到的TWSC结果中出现了显著的南北条带误差，所以必须对GRACE时变重力场模型进行滤波处理以削弱上述误差.常用的滤波算法有高斯滤波[6]、Fan滤波[7]、去相关滤波[8-9]、Han滤波[10]等，不同滤波算法的处理效果不尽相同，并直接影响到反演结果.稀慧等[11]采用多种滤波方法推算全球平均海水质量变化，与去相关滑动滤波的比较结果显示，由不同滤波算法处理得到全球平均海水质量变化结果具有显著的区别；张青全等[12]比较了4种不同滤波算法在GRACE反演西南地区岩溶区陆地水储量变化的影响，指出不同滤波反演结果在空间分布上差异较大.但是上述比较是以其他滤波方法结果和降水数据为基准，并不能反映出真实的陆地水储量变化.

本研究采用2002年4月至2017年6月连续15年的GRACE时变重力场模型数据反演全球陆地水储量变化，并利用高斯滤波、Fan滤波、去相关滤波和Han滤波分别对相关误差进行处理.为了对比4种滤波算法的处理结果，本研究从时空分布，长期趋势变化和季节性变化3个方面分别进行阐述，并将4种滤波算法结果与NASA提供的结果进行比较.

1 数据与算法

1.1 GRACE数据处理

采用由美国德克萨斯大学空间研究中心提供的2002年4月至2017年6月的GRACE RL06时变重力场球谐系数，其截断阶数为60.在进行反演前需要对模型的球谐系数进行一系列的预处理，具体步骤如下：1)采用卫星激光测距(Satellite Laser Ranging,SLR)获取的高精度C20项数据对重力场模型相应项的系数进行替换[13]；2)利用文献[14]的成果对模型一阶项进行地心变化改正；3)采用滤波算法对相关误差进行处理.

利用GRACE时变重力场反演TWSC的计算公式如下，用等效水高(Equivalent Water High,EWH)表示TWSC[7]，有，

ΔSlmsin(mλ))

(1)

1.2 滤波算法原理

1.2.1 高斯滤波

该滤波是将每个点的密度变化用所有点密度变化的加权平均值替代，以达到平滑效果.当式(1)中Wlm=Wl时，即空间平滑函数只与阶数相关，可根据递归关系计算高斯滤波系数Wi(i=1,…,l)，具体如下：

(2)

式中，b=ln2/[1-cos(r/a)]；r为滤波半径/km;e为自然常数.

1.2.2 Fan滤波

该算法本质上是各向异性滤波，即对球谐系数的阶数和次数均采用与高斯滤波相同的处理方式，即式(1)中Wlm=Wl·Wm，具体表达式如下:

ΔSlmsin(mλ))

(3)

式中，Wj(j=1,…,m)的求解方法与式(2)完全一致.

1.2.3 各向异性高斯滤波

考虑到GRACE时变重力场模型误差是各向异性[10]，因此，高斯滤波仅对阶相关误差进行处理的方法存在着缺陷，各向异性高斯滤波系数的具体表达式如下：

(4)

式中，r0和r1为滤波半径,/km，且r0

1.2.4 去相关滤波

去相关滤波的基本原理是保持重力场模型前l×l阶的位系数不变，将n阶多项式拟合得到的大于或等于m阶次位系数从原模型中扣除，以消除误差[15-16].

1.3 线性拟合模型

为了详细分析TWSC在时域中的变化规律，一般采用线性自回归方程从TWSC时间序列中提取长期趋势项、长期趋势加速度项、周年项和半周年项，其具体表达式如下[17]：

ΔEWH(t)=a+bt+ccos(2πt)+dsin(2πt)+ecos(4πt)+fsin(4πt)+ε

(5)

式中，a为常数，b为长期趋势，c和d为周年，e和f为半周年，ε为残差，t为时间.采用最小二乘配置法求解上述方程.同时，TWSC时间序列中的周年项和半周年项的振幅(Aann和Asemi-ann)、相位(Φann和Φsemi-ann)计算公式如下[18]：

(6)

2 实验结果与分析

本研究计算了2014年10月全球陆地水储量变化，结果如图1所示.由图1(a)可明显看出陆地水储量变化结果呈现出显著的南北分布条带误差，由于误差的干扰造成了正常信号提取的困难，因此对条带误差的处理是必要的.分别采用350 km 高斯滤波，300 km Fan滤波，各向异性高斯滤波和P3M6多项式滤波对条带误差进行了处理，结果如图1(b)～(e)所示.对比可知，Fan滤波和各向异性高斯滤波的处理效果最为显著，几乎看不到条带误差的痕迹，但是在削弱条带误差影响的同时，真实信号也受到了影响.比较图(c)和(d)可知，在相同误差处理效果的前提下，各向异性高斯滤波比Fan滤波能更多地保留真实信号.结合图(b)和(e)可知，P3M6多项式滤波处理中低纬度地区的条带误差处理效果较好，而高斯滤波则对高纬度地区的误差处理效果较好.综上所述，各向异性高斯滤波既能很好地处理条带误差的影响，又可以最大限度的保留真实信号.

图1 全球陆地水储量变化滤波结果(2014年10月)

为了进一步比较上述4种滤波算法的处理效果，本研究计算了2002年4月至2017年6月的全球陆地水储量变化的时间序列如图2(a)～(d)所示.由于4种算法的结果非常接近，由图可知4组时间序列的变化趋势基本一致，同时在数量级上，除了Fan滤波结果以外，其他滤波也完全相同.比较图2(b)和图2(a)、(c)、(d)发现，Fan滤波结果要比其他3种滤波算法结果小一个数据量级，这再次印证了Fan滤波除了削弱条带误差之外，对真实信号的削弱效果也较其他3种滤波算法更为显著.本研究计算了4组时间序列的长期趋势变化、周年振幅、周年相位、半周年振幅和半周年相位，结果如表1所示.

图2 全球陆地水储量变化时间序列(2002年4月至2017年6月)

由表1可知，在长期趋势变化、周年振幅、周年相位、半周年振幅和半周年相位上，350 km 高斯滤波、300 km Fan滤波和P3M6多项式滤波3者结果较为接近.而各向异性高斯滤波则在长期趋势变化、周年振幅和半周年振幅三方面存在着一定的差别，但是考虑到3种指标的数量级都很小，这种差别是可以忽略不计的.在周年相位和半周年相位两种指标方面，350 km 高斯滤波、300 km Fan滤波和各向异性高斯滤波的结果一致，但是300 km Fan滤波与上述3种滤波算法的结果则存在着一定的差异，这与图2所得的结果相互印证.

表1 全球陆地水储量变化长期趋势变化、周年项和半周年项

3 结 论

针对GRACE时变重力场模型反演陆地水储量变化中出现的条带误差问题，本研究采用4种常见的滤波算法(高斯滤波、Fan滤波、各向异性高斯滤波和P3M6多项式滤波)分别对条带误差进行处理.为了比较4种滤波算法的误差处理效果，分别比较了4种滤波算法计算得到的全球陆地水储量变化及其时间序列，以及长期变化趋势、周年变化项和半周年变化项.结果表明，各向异性高斯滤波和Fan滤波在进行误差处理方面要优于其他两种算法，在保留真实信号方面，各向异性高斯滤波要优于Fan滤波;在长期变化趋势、周年变化和半周年变化方面，除了Fan滤波以外，其他3种滤波算法的结果基本上完全一致.综上所述，各向异性高斯滤波更适合用于对GRACE模型条带误差的处理.