GRACE滤波数据处理之DDK系列滤波

以CSR RL06无约束解为例，进行DDK1-8滤波数据处理，人为构造如下有关读取数据的控制文件：

Github上下载DDK滤波核函数：GitHub - strawpants/GRACE-filter: Contains software for filtering (destriping) GRACE Stokes coefficientshttps://github.com/strawpants/GRACE-filter

close all;clearvars -except;
addpath('E:\1A\24\Code\Tool\Function');
addpath('E:\1A\24\Code\Tool\Function\GRACE_functions');
fid=fopen('E:\1A\24\Code\txt\Contr_CSR.txt','r');
num_file= fscanf(fid,'%d',1);
dir_c20 = fscanf(fid,'%s',1);dir_degree_1= fscanf(fid,'%s',1);
dir_cs21= fscanf(fid,'%s',1);dir_cs22= fscanf(fid,'%s',1);
dir_in = fscanf(fid,'%s',1);dir_out= fscanf(fid,'%s',1);
FileNameTime02={'2002-04','2002-05','2002-08','2002-09','2002-10','2002-11','2002-12'};
FileNameTime03={'2003-01','2003-02','2003-03','2003-04','2003-05','2003-07','2003-08','2003-09','2003-10','2003-11','2003-12'};
FileNameTime04={'2004-01','2004-02','2004-03','2004-04','2004-05','2004-06','2004-07','2004-08','2004-09','2004-10','2004-11','2004-12'};
FileNameTime05={'2005-01','2005-02','2005-03','2005-04','2005-05','2005-06','2005-07','2005-08','2005-09','2005-10','2005-11','2005-12'};
FileNameTime06={'2006-01','2006-02','2006-03','2006-04','2006-05','2006-06','2006-07','2006-08','2006-09','2006-10','2006-11','2006-12'};
FileNameTime07={'2007-01','2007-02','2007-03','2007-04','2007-05','2007-06','2007-07','2007-08','2007-09','2007-10','2007-11','2007-12'};
FileNameTime08={'2008-01','2008-02','2008-03','2008-04','2008-05','2008-06','2008-07','2008-08','2008-09','2008-10','2008-11','2008-12'};
FileNameTime09={'2009-01','2009-02','2009-03','2009-04','2009-05','2009-06','2009-07','2009-08','2009-09','2009-10','2009-11','2009-12'};
FileNameTime=[FileNameTime02,FileNameTime03,FileNameTime04,FileNameTime05,FileNameTime06,FileNameTime07,FileNameTime08,FileNameTime09];

FileNameTimeChar=char(FileNameTime);
int_year=zeros(num_file,1);int_month=zeros(num_file,1);%用于存储年和月
for i=1:num_file
    int_year(i)=str2double(FileNameTimeChar(i,1:4));
    int_month(i)=str2double(FileNameTimeChar(i,6:7));
end
time=int_year+(int_month-0.5)/12;
% save E:\1A\24\Result\time.mat time int_year int_month FileNameTime
file_name=cell(num_file,1);
for i = 1:num_file
    file_name{i} = fscanf(fid,'%s',1);
end
fclose(fid);
name=file_name{1,1}(1,30:32);%CSR用这个
lmax=60;
cs= zeros(num_file,lmax+1,lmax+1);
cs_sgi= zeros(num_file,lmax+1,lmax+1);
cs_res= zeros(num_file,lmax+1,lmax+1);
cs_mss= zeros(num_file,lmax+1,lmax+1);

tic;
hwait=waitbar(0,'Waiting>>>>>>>>');  %加载等待对话框
 for ii=1:num_file
     str=['Processing...',num2str(ii),'/',num2str(num_file),'    '];
     hwait=waitbar(ii/num_file,hwait,str,'Name','SSM');
 pathname=strcat(dir_in,file_name{ii,1});
 [cs(ii,:,:),cs_sgi(ii,:,:)] =gmt_readgfc(pathname);
 end
toc

cs_replace=cs;
[cs_replace] = gmt_replace_degree_1(dir_degree_1,cs_replace,int_year,int_month,num_file);
[cs_replace] = gmt_replace_C20(dir_c20,cs_replace,int_year,int_month,num_file);
[cs_replace] = gmt_replace_C30(cs_replace,int_year,int_month,num_file);

cs_mean = mean(cs_replace(19:90,:,:)); %求2004年1月至2009年12月间的平均值

for i=1:num_file
    cs_res(i,:,:)  = cs_replace(i,:,:)-cs_mean(1,:,:);%扣除平均值
end

%%DDK滤波
type_filter = 'DDK3';
addpath('E:\1A\24\Code\Tool\ddk_filtercoef');
switch type_filter
    case 'DDK1'
        path = 'Wbd_2-120.a_1d14p_4';
    case 'DDK2'
        path = 'Wbd_2-120.a_1d13p_4';
    case 'DDK3'
        path = 'Wbd_2-120.a_1d12p_4';
    case 'DDK4'
        path = 'Wbd_2-120.a_5d11p_4';
    case 'DDK5'
        path = 'Wbd_2-120.a_1d11p_4';
    case 'DDK6'
        path = 'Wbd_2-120.a_5d10p_4';
    case 'DDK7'
        path = 'Wbd_2-120.a_1d10p_4';
    case 'DDK8'
        path = 'Wbd_2-120.a_5d9p_4';
end
% 开始滤波
W = read_BIN(path);
for ii = 1:num_file
    sc = gmt_cs2sc(reshape(cs_res(ii,:,:),[lmax+1 lmax+1]));
    s = sc(:,1:lmax+1);
    s(:,end) = 0;
    s = fliplr(s);
    c = sc(:,61:end);
    [c_filt,s_filt] = filterSH(W,c,s);%此处有滤波函数
    s_filt = fliplr(s_filt);
    sc_filt = [s_filt(:,1:end-1),c_filt];
    cs_filt = gmt_sc2cs(sc_filt);
    cs_res(ii,:,:) = cs_filt;
end

radius_filter=0;
for i=1:num_file
cs_tmp1(:,:) = cs_res(i,:,:);
cs_tmp2(:,:)=gmt_gc2mc(cs_tmp1);      %%球谐系数转换为质量系数 
cs_tmp3(:,:) = gmt_destriping(cs_tmp2,'NONE');%CHENP4M6 NONE  SWENSON 去相关滤波
cs_mss(i,:,:)=gmt_gaussian_filter(cs_tmp3,radius_filter);
end
% eval(['cs_mss' type_filter '=cs_mss;']);

eval(['grid_' name '=gmt_cs2grid(cs_mss,0,1);']);
% save cs_mssDDK.mat cs_mssDDK1 cs_mssDDK2 cs_mssDDK3 cs_mssDDK4 ...
%     cs_mssDDK5 cs_mssDDK6 cs_mssDDK7 cs_mssDDK8 

ii=1;
eval(['grid_1(:,:)=grid_' name '(:,:,ii);']);
ii=1;
eval(['grid_1(:,:)=grid_' name '(:,:,ii);']);
gmt_grid2map(grid_1*100,-30,30,1,0,'cm',[name   ' ' num2str(int_year(ii),'%02d') '-' num2str(int_month(ii),'%02d') ' ' type_filter],20)

参考文献：

1、Kusche J. Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models[J]. Journal of Geodesy, 2007, 81(11): 733-749.

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

欢迎多多交流