论文推荐｜郭飞霄：GRACE时变重力场各向异性高斯组合滤波方法

《测绘学 》

构建与学术的桥梁 拉近与权威的距离

GRACE时变重力场各向异性高斯组合滤波方法

郭飞霄^1,2,3

1. 信息工程大学, 河南 郑州 450001;

2. 地理信息工程国家重点实验室, 陕西 西安 710054;

3. 西安测绘研究所, 陕西 西安 710054;

4. 西安航天天绘数据技术有限公司, 陕西 西安 710054;

5. 31002部队, 北京 100094

收稿日期：2018-08-22；修回日期：2019-01-25

基金项目：国家自然科学基金（41674082；41774018）；大地测量与地球动力学国家重点实验室开放基金（SKLGED2017-3-2-E）

关键词：组合滤波 条带噪声 各向异性高斯滤波 均方根滤波 信噪比

Non-isotropic Gaussian combination filtering method of GRACE time-variable gravity

GUO Feixiao^1,2,3, SUN Zhongmiao^2,3, REN Feilong⁴, WANG Feifei⁵

1. Information and Engineering University, Zhengzhou 450001, China;

2. State Key Laboratory of Geo-Information Engineering, Xi‘an 710054, China;

3. Xi‘an Research Institute of Surveying and Mapping, Xi‘an 710054, China;

4. Xi‘an Aerospace Data Technology Co., Ltd, Xi‘an 710054, China;

5. Troops 31002, Beijing 100094, China

Foundation support: The National Natural Science Foundation of China (Nos. 41674082;41774018); The Open Foundation of State Key Laboratory of Geodesy and Earth’s Dynamics (No. SKLGED2017-3-2-E)

First author: GUO Feixiao(1988—), male, PhD candidate, engineer, majors in satellite gravimetry data processing and application.E-mail: gravity_gfx@126.com.

Abstract: Affected by factors of measurement errors and so on, the surface mass variation inversion results would be dominated by serious north-to-south stripe noise with using time-variable gravity models directly. Filtering methods must be adopted. Optimal parameters of different filtering methods are determined by taking signal-to-noise ratio maximum as the criterion. On this basis, a new non-isotropic combination filtering method is proposed based on the features of time-variable gravity model spherical harmonic serrors. The new non-isotropic combination filtering method is different from traditional two-steps combination filtering methods. It just performs in one step with combining non-isotropic Gaussian filtering and RMS filtering. The basic idea of the new combination filtering method is to apply larger weights to low-order spherical harmonic coefficients while smaller weights to high-order spherical harmonic coefficients for keeping more signals and suppressing noise. The experiment results show that the new combination filtering method is simple to calculate and stripe noise can be removed effectively. The new combination filtering method improves the signal-to-noise ratio compared with previous filtering methods and traditional two-steps combination filtering methods.

Key words: combination filtering stripe noise non-isotropic gaussian filtering root mean square filtering signal-to-noise ratio

随着科学技术不断发展，人类对地球认识水平不断提高，有效监测地球环境变化是当今地球科学研究的重要任务之一。地球环境变化与地表物质迁移密切相关，全球海平面上升、冰川消融、陆地水储量变化等全球性环境变化问题，都与地表物质迁移过程有着密切联系^[1-2]。因此，开展地表物质迁移研究对监测地球环境变化等具有重要意义。GRACE(gravity recovery and climate experiment)卫星于2002年3月成功发射，开创了高精度全球重力场观测与气候变化试验的新纪元，为连续监测地表物质迁移和重新分布提供了直接观测手段，目前已广泛应用于陆地水储量变化、冰川融化和海水质量变化反演等领域^[3-5]。

球谐系数法是利用卫星时变重力场模型反演地表物质迁移研究中应用最广泛的方法。但是，由于卫星载荷的仪器测量误差、混频误差及卫星轨道等因素的影响^[6]，直接使用时变重力场模型的球谐系数法反演结果存在严重南北条带噪声，必须对球谐系数进行滤波处理。按照滤波思想的不同主要有两类方法，一类方法是高斯滤波^[7]、扇形滤波^[8]、维纳滤波^[9]、均方根(root mean square，RMS)滤波^[10]等，这类方法基本思想是通过降低模型高阶次项系数权重以减小短波长分量误差影响；另一类方法则是通过拟合多项式消除球谐系数奇偶次项间的相关误差，主要有滑动多项式拟合方法^[11]、PnMm方法^[12]及反向延拓滑动多项式拟合方法等^[13]。常用的高斯滤波和扇形滤波方法在低纬度地区具有较好的去条带滤波效果，而去相关滤波则在高纬度地区具有较好的去条带效果^[4]。因此，为了最大限度消除条带噪声并充分发挥不同滤波方法的各自优势，国外学者提出了组合滤波法，最普遍的组合方式是将高斯滤波或扇形滤波等方法与去相关滤波进行两步滤波处理，其基本思想为：首先对球谐系数进行去相关滤波，分别消除球谐系数奇数和偶数次项之间的相关误差；再对经去相关滤波后的系数进行高斯滤波等滤波处理，削弱高阶次项噪声。研究表明：组合滤波同时削弱了球谐系数高阶次项误差和奇偶次项间相关误差，相比单一滤波方法，去条带噪声效果更好^[14]。此外，还有学者提出了两步法各向异性组合滤波^[15]，该方法分两步分别进行各向异性高斯滤波和RMS滤波，保留了低阶项真实信号并压制高阶项噪声，滤波效果较好。但是，当前的组合滤波方法实质上都是两步处理，不同步骤的滤波方法选取以及滤波参数选取都具有一定经验性，处理较为复杂。

1 滤波方法比较分析

文献[7]最早提出了球谐系数法，该方法将研究时间段内所有时变重力场模型的平均值作为稳态重力场模型，通过计算相对于稳态重力场模型球谐系数的改变量进而反演得到地表物质迁移，即

(1)

式中，r_a为地球平均半径；θ、λ为地心余纬和地心经度；l、m为球谐函数的展开阶数和次数；ΔClm、ΔSlm为球谐系数改变量；Plm(cosθ)为完全规格化的勒让德缔合函数；kl为负荷勒夫数；W_l为滤波系数；ρ_w为水密度；ρ_a为地球平均密度；Δh为质量变化的等效水柱高。

1.1 高斯滤波

文献[19]最早提出空间平滑函数，其本质是将某点的密度变化看成所有点密度变化的加权平均以达到平滑效果。文献[7]首先采用归一化高斯平均核函数，建立了高斯滤波方法，并根据递归关系来计算高斯滤波系数

(2)

式中，r是滤波半径。

1.2 各向异性高斯滤波

研究表明：GRACE时变重力场模型高阶高次项系数误差大于高阶低次项系数误差，即时变重力场模型误差是各向异性的^[20-21]。因此，仅采用阶相关的高斯滤波会压制高阶低次项的真实信号。针对此，文献[20]提出了滤波系数与所用模型数据的阶和次数都相关的各向异性高斯滤波

(3)

式中，W为高斯滤波核函数；r₁和r₀为滤波半径且r₀小于r₁；m₁为选定的次，且0≤m≤m₁。

1.3 扇形滤波

文献[8]提出了扇形滤波，该方法对球谐系数的阶和次项同时进行高斯滤波处理

(4)

式中，W_l和W_m分别为与阶、次数相关的高斯滤波系数。由于其滤波系数在阶-次平面上投影形状为扇形，故称扇形滤波。

1.4 维纳滤波

高斯滤波、各向异性高斯滤波及扇形滤波都需要选取合适的滤波半径，滤波半径的选取直接决定了反演结果中条带噪声的去除效果。滤波半径越大，去条带噪声效果越好，但对模型空间分辨率牺牲程度也越大；滤波半径过小，则无法有效消除条带噪声，给真实信号识别带来困难。文献[9]基于信号和噪声的阶方差谱来建立信号和噪声函数，建立了维纳滤波

(5)

式中，σs,l²和σn,l²分别是信号和噪声的阶方差。该方法通过设计滤波函数，对测量信号进行线性卷积得到实际输出信号，使其与期望输出信号满足最小二乘，进而求解得到滤波系数，其滤波系数只依赖于信号和噪声模型的阶方差，属于各向同性滤波。

1.5 RMS滤波

由于GRACE卫星观测到的时变重力变化信号在陆地明显高于在海洋，根据该条件可以陆地和海洋信噪比最大为准则建立最优滤波。文献[10]基于该准则得到滤波系数，建立了RMS滤波

(6)

式中，MDL为水文模型加上海洋模型的球谐系数，GRC为GRACE时变重力场模型球谐系数，a、b为滤波参数，MASS为质量变化，Err为GRACE测量误差，其中a、b选取满足信噪比RMS_Ratio达到最大值

(7)

1.6 两步法组合滤波

由于去相关滤波仅对高纬度地区具有较好去条带效果，低纬度地区条带噪声无法消除^[4]。因此，去相关滤波不能单独使用，需要与其他滤波方法组合使用。当前，最普遍的组合滤波方式是将高斯滤波或扇形滤波与去相关滤波分两步进行滤波处理，即首先采用去相关滤波对球谐系数奇数和偶数阶项之间的相关误差分别进行消除，再对经过去相关滤波后的球谐系数进行高斯滤波或扇形滤波。

1.7 滤波比较分析

试验采用GFZ Release 05版本时变重力场模型数据，时间跨度为2004年1月至2010年12月。反演时球谐系数截断至60阶，C₂₀项采用激光测距数据所得结果替换，RMS滤波中的模型数据分别选GLDAS和ECCO模型。试验对不同滤波参数条件下的滤波结果信噪比值进行了统计，以信噪比最大为准则，对滤波参数的最优选取进行了分析。由于不同月份滤波结果类似，受篇幅所限，下文以2010年1月全球陆地水储量变化结果为例进行讨论分析。

图 1所示分别为高斯滤波和扇形滤波结果信噪比值随滤波半径变化。从图中结果可知，反演结果信噪比值均大于1，信噪比值随滤波半径增大先增大后减小；当滤波半径为500 km时，高斯滤波结果信噪比值达到最大值；滤波半径为400 km时，扇形滤波结果信噪比值达到最大值。因此，后续试验中高斯滤波半径取500 km，扇形滤波半径取400 km。

图 1 不同滤波半径滤波结果信噪比值Fig. 1 Ratio of filtering results in different radius

图选项

当时变重力场模型球谐系数的次数大于等于15时，球谐系数误差越来越显著^[13]，故试验中各向异性高斯滤波计算公式(3)中取m₁=15。根据图 1结果，滤波半径为500 km时，高斯滤波结果信噪比值达到最大，结合1.2节各向异性高斯滤波系数公式，试验中各向异性高斯滤波半径分别取r₀=200、300、400和500 km及r₀<r₁≤1000 km。表 1所示为不同滤波半径组合情况下，各向异性高斯滤波结果信噪比值统计，从表 1结果可看出，当r₀=200 km、r₁=500 km时，各向异性高斯滤波结果信噪比达到最大值。因此，后续试验中各向异性高斯滤波半径r₀=200 km、r₁=500 km。

表 1 不同滤波半径组合各向异性高斯滤波结果信噪比统计Tab. 1 Ratio statistics of non-isotropic Gaussian filtering in different radius combinations

r0	r1
	300	400	500	600	700	800	900	1000
200	2.827	2.933	2.934	2.915	2.889	2.860	2.831	2.801
300	2.863	2.881	2.865	2.840	2.811	2.781	2.752
400	2.796	2.786	2.761	2.733	2.703	2.674
500	2.715	2.692	2.663	2.634	2.605

表选项

图 2所示为RMS滤波参数a、b分别取1至3时滤波结果信噪比值。从图 2结果可看出，当a=2、b=1.8时，信噪比值达到最大。因此，后续试验中RMS滤波最优参数取为a=2、b=1.8。

图 2 不同a、b取值滤波结果信噪比值Fig. 2 Ratio of filtering results in differentaandb

图选项

记组合滤波1为高斯滤波(r=300 km)+去相关滤波，组合滤波2为高斯滤波(r=400 km)+去相关滤波，组合滤波3为高斯滤波(r=500 km)+去相关滤波，组合滤波4为扇形滤波(r=300 km)+去相关滤波，组合滤波5为扇形滤波(r=400 km)+去相关滤波，组合滤波6为扇形滤波(r=500 km)+去相关滤波。

图 3为不同滤波方法所得全球陆地水储量变化反演结果。从图中可以看出，仅采用去相关滤波的反演结果条带噪声严重，其他滤波方法则对条带噪声均具有较好消除效果，各向异性高斯滤波和RMS滤波所得结果信号较强，组合滤波1结果在低纬度海洋区域上存在较为明显条带噪声，组合滤波2结果也在低纬度海洋区域残存少量条带噪声。

图 3 不同滤波方法比较Fig. 3 Comparison of different filtering methods

图选项

表 2为不同滤波方法结果信噪比值统计。从表 2结果可以看出，各向异性滤波结果信噪比值最大，其次分别是组合滤波4、RMS滤波，而组合滤波1结果信噪比值最小；并且，维纳滤波结果信噪比略大于高斯滤波(r=500 km)结果信噪比，与扇形滤波(r=400 km)结果信噪比相当。以上结果表明：单一滤波方法中，各向异性高斯滤波和RMS滤波结果信噪比值最大，在消除条带噪声的同时保留了更多有效真实信号；并且，组合滤波结果的信噪比值大于单一滤波方法结果。

表 2 不同滤波方法结果信噪比值Tab. 2 Ratio of filtering results with different methods

滤波方法	信噪比值
高斯滤波(r=500 km)	2.641
各向异性高斯滤波(r0=200 km、r1=500 km)	2.934
扇形滤波(r=400 km)	2.754
维纳滤波	2.762
RMS滤波(a=2、b=1.8)	2.778
组合滤波1	2.359
组合滤波2	2.695
组合滤波3	2.683
组合滤波4	2.836
组合滤波5	2.769
组合滤波6	2.670

表选项

2 各向异性组合滤波

由上文试验结果可知，单一滤波方法中，各向异性高斯滤波和RMS滤波结果信噪比值最大。图 4(a)和(b)所示分别为各向异性高斯滤波(r₀=200 km、r₁=500 km)和RMS滤波(a=2、b=1.8)各阶次滤波系数值。从图 4可以看出，各向异性高斯滤波系数在阶数确定的条件下随着次数的增大而减小，对精度较差的高阶高次项权重较小，可有效压制高阶次项噪声；与各向异性高斯滤波相比，RMS滤波则对精度较高的低次项权重更大，更好保留了低次项信号，并且对高阶高次项滤波系数权重也较小，也可有效压制高阶次项噪声。

图 4 不同滤波方法C_nm项各阶次滤波系数值Fig. 4 Filtering coefficient values of different filtering methods

图选项

(8)

式中，Wlm¹和Wlm²分别为各向异性高斯滤波和RMS滤波系数值。

3 各向异性组合滤波试验分析

试验中时变重力场模型数据和预处理过程同1.7节。各向异性组合滤波分界次数m_k分别取5、10、15和20，以2010年1月全球陆地水储量变化反演结果为例，所得结果如图 5所示。从图 5可以看出，4种情形下的各向异性组合滤波整体上均有效消除了条带噪声，但随着分界次数m_k增大，各向异性组合滤波结果在海洋区域上逐渐呈现微弱条带噪声。

图 5 不同参数各向异性组合滤波结果Fig. 5 Ground water storage variation inversion results using non-isotropic combination filtering in different parameters

图选项

表 3 不同参数各向异性组合滤波信噪比值Tab. 3 Ratio of non-isotropic combination filtering in different parameters

各向异性组合滤波	信噪比值
mk=5	2.949
mk=10	2.958
mk=15	2.940
mk=20	2.920

表选项

图 6 采用不同阶模型的各向异性组合滤波结果Fig. 6 Ground water storage variation inversion results using non-isotropic combination filtering in different degree of models

图选项

表 4 不同阶数模型组合滤波结果信噪比值Tab. 4 Ratio of non-isotropic combination filtering results in different degrees

模型阶数	信噪比值
20	2.757
30	2.860
40	2.921
50	2.947
60	2.940

表选项

4 结论

(1) 在信噪比最大准则条件下，高斯滤波最优半径为500 km，各向异性高斯滤波最优半径为r₀=200 km、r₁=500 km，扇形滤波最优半径为400 km，RMS滤波最优参数为a=2、b=1.8；维纳滤波结果略优于高斯滤波(r=500 km)结果，与扇形滤波(r=400 km)去条带效果相当；单一滤波方法中，各向异性高斯滤波(r₀=200 km、r₁=500 km)和RMS滤波(a=2、b=1.8)的反演结果信噪比值最大。

由于国内外对GRACE卫星时变重力场模型数据滤波效果的评判尚无公认准则，仅采用信噪比评判滤波效果还不够全面。而滤波的同时往往也会造成信号泄露等问题，因此滤波方法选择和滤波参数选取应当在真实信号衰减和噪声减小之间作出平衡。后续，关于条带噪声产生机理、球谐系数误差特性、信噪比与泄露误差间关系、滤波效果评判准则等问题还需作进一步的研究。

【引文格式】郭飞霄, 孙中苗, 任飞龙, 等. GRACE时变重力场各向异性高斯组合滤波方法[J]. 测绘学 ，2019，48(7)：898-907. DOI: 10.11947/j.AGCS.2019.20180365

权威 | 专业 | 学术 | 前沿

微信投稿邮箱 | song_qi_fan@163.com

进群请备注：姓名+单位+稿件编号