GRACE数据的滤波方法实现（Han、Fan和Wiener）

本文基于GRACE_Matlab_Toolbox工具包的基础上利用matlab实现一些其他的滤波函数，例如非各向同性滤波（Han）、扇形滤波（Fan）和维纳滤波（Wiener），在介绍上述滤波之前，先分析下高斯滤波：

高斯滤波系数递归公式如下：

GRACE_Matlab_Toolbox工具包利用近似公式计算高斯滤波系数，公式和代码如下：

W = gmt_gaussian(lmax,radius_filter);

1、Han滤波

clear;
max_degree=60;
W = zeros(max_degree+1,1);
r0=200;r1=300;m1=15;
for i=1:max_degree+1
    for j=1:i
      rm=r0+(r1-r0)*j/m1;
      W(i,j)=exp(-((i-1)*rm/6371)^2/(4*log(2)));
    end
end

%将非各向同性滤波系数扩展为61*121
w=zeros(61,121);lmax=60;
w(:,61:121)=W;jj=2;
for ii=lmax:-1:1
    w(:,ii)=W(:,jj);jj=jj+1;
end

m1=0:60;l1=0:60;
[m,l]=meshgrid(m1,l1);
figure(1)
pcolor(m,l,W);shading flat;
caxis([0,1]);colormap('jet');grid off;
xlabel('Order(m)','Fontname','Time New Roman',"FontSize",20);
ylabel('Degree(l)','Fontname','Time New Roman',"FontSize",20);

示例代码中滤波半径r0=200km、r1=300km、选定的次m1=15，绘制滤波系数图如下

2、Fan滤波

扇形滤波是对阶数和次数同时进行高斯滤波：

针对GRACE_Matlab_Toolbox工具包的gmt_gaussian_filter函数进行修改如下：

Wl = gmt_gaussian(lmax,radius_filter);%阶数高斯滤波系数
Wm = gmt_gaussian(lmax,radius_filter);%次数高斯滤波系数
sc(:,:) = gmt_cs2sc(field);% fan filter(gauss filter)
%先对对次数进行高斯滤波
for ml = lmax:2*lmax
    sc(:,ml+1) = Wm(ml-lmax+1)*sc(:,ml+1);
end
for ml = lmax:-1:0
    sc(:,ml+1) = Wm(lmax-ml+1)*sc(:,ml+1);
end
%再对阶次数进行高斯滤波
for ll = 0:1:lmax
    scnew(ll+1,:) = Wl(ll+1)*sc(ll+1,:);
end

3、Wiener滤波

主要依靠三个公式计算滤波系数h：

首先以等效水柱高形式绘制测量信号的阶方差图像，再运用上述前两个公式进行最小二乘拟合，求得拟合参数a、b、c、d，最后利用第三公式计算滤波系数h（注意：对于最大阶数60的球谐产品，上述第一个公式适用于0-20阶，第二个公式适用于31-60阶数）

%0-20阶
l=1:60;
Z1(2:20,1)=log10(Temp(2:20,18));L1(1:20,1)=log10(l(1,1:20))';
% 
p = polyfit(L1(2:20,1), Z1(2:20,1), 1);b=-p(1,1);a=p(1,2);
for j=1:lmax
    y11(j,1)=10.^a/(j.^b);
end
Z11=-b*L1(2:20,1)+a;

%30阶以后
Z2(:,1)=log10(Temp(31:60,18));L2(:,1)=l(1,31:60)';
p1 = polyfit(L2, Z2, 1);d=p1(1,1);c=p1(1,2);
for j=1:lmax
    y21(j,1)=10.^(c+d*j);
end

figure(1);
l=2:60;th=1.5;
plot(l,(Temp(2:end,18)),'color',[0.2 0.6 0.2],'Linewidth',th);  
set(gca,'yscale','log');
xticks([0 10 20 30 40 50 60]);
xlabel('Degree','Fontname','Time New Roman');
ylabel('Geoid degree error(m)','Fontname','Time New Roman');
grid on; hold on;

plot((1:60)',y11,'color',[1 0.6 0.2],'Linewidth',th);
plot((1:60)',y21,'color',[0 1 0],'Linewidth',th);
h=y11./(y11+y21);plot(1:60,h'); %wiener滤波系数

两段拟合情况如下：

绘制2003年11月CSR RL06 Gaussian、Han、Fan、Wiener滤波结果： 

参考文献：

[1] 郭飞霄. 地表物质迁移的卫星大地测量反演理论与方法研究[D].战略支援部队信息工程大学,2019.DOI:10.27188/d.cnki.gzjxu.2019.000012.

[2] Han S C, Shum C K, Jekeli C, et al. Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement[J]. Geophysical Journal International, 2005, 163(1): 18-25.

[3] Sasgen I, Martinec Z, Fleming K. Wiener optimal filtering of GRACE data[J]. Studia Geophysica et Geodaetica, 2006, 50: 499-508.

[4] Chambers D P. Observing seasonal steric sea level variations with GRACE and satellite altimetry[J]. Journal of Geophysical Research: Oceans, 2006, 111(C3).