本人:B站同步

我们知道，任何球面上的一组观测值可以球谐展开至一定的阶数，见如下的公式：

而为了将GLDAS提取的水储量与GRACE能够进行对比，需要进行同样的截断滤波处理。本专栏以GLDAS数据的土壤水为例，分析如何进行球谐展开并滤波。主要代码如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLDAS expand 
% Chistrong Wen
% 2022-12-27
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GLDAS INFO
file = 'GLDAS_VIC10_M.A201710.021.nc4';
ncdisp(file)

%% READ SOIL MOISTURE
lon = ncread(file,'lon');
lat = ncread(file,'lat');
[lon,lat] = meshgrid(lon,lat);
%% 
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);

%% SoilMoi0_30cm_inst SoilMoi_depth2_inst SoilMoi_depth3_inst
soil1 = ncread(file,'SoilMoi0_30cm_inst');
soil2 = ncread(file,'SoilMoi_depth2_inst');
soil3 = ncread(file,'SoilMoi_depth3_inst');
soil_sum = soil1+soil2+soil3;
out_soil.lon = lon;
out_soil.lat = lat;
fill = nan(30,360);
out_soil_temp = [fill;soil_sum'];
out_soil.rg  = out_soil_temp;
figure
wzq_plot(out_soil)
title('raw data')
grid = out_soil.rg;

degree_max = 60;
[cs] = gmt_grid2cs(grid,degree_max);

destrip_method = 'SWENSON';
radius_filter = 300; 
type = 1;

[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);

figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg  = grid_filter;
wzq_plot(out_soil_filter)
title('degree 60')

%% truncted to 96
degree_max = 96;
[cs] = gmt_grid2cs(grid,degree_max);

destrip_method = 'SWENSON';
radius_filter = 300; 
type = 1;

[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);

figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg  = grid_filter;
wzq_plot(out_soil_filter)
title('degree 96')
%% smooth at 500 km
degree_max = 60;
[cs] = gmt_grid2cs(grid,degree_max);

destrip_method = 'SWENSON';
radius_filter = 500; 
type = 1;

[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);

figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg  = grid_filter;
wzq_plot(out_soil_filter)
title('Gaussian 500 km')

下面是截断滤波后的结果：我们看到截断滤波后，信号强度有所衰减，而随着截断阶数越大，信号越接近真实信号强度。而高斯平滑的半径越大，信号也衰减。

欢迎大家收藏点赞❤❤❤