球谐系数阶方差计算和图谱绘制

运用下面这期博文中的格网数据

绘制全球各大洲典型流域的时间序列图-CSDN博客https://blog.csdn.net/WH_G_Y/article/details/139972124?spm=1001.2014.3001.5501大地水准面阶方差计算公式：

 计算中所用到两个转换函数：

1、“GenerateC2030C2021TransIndex”函数  功能：调换球谐系数的排列顺序（例如按阶或者按次升序排列）

function [C2021toC2030_index,C2030toC2021_index]=GenerateC2030C2021TransIndex(lmax)
% The CSR SHC storage: C20,C30,C40,...,C21,S21,C22,S22,C31,S31,C32,S32,C33,S33,...
% We define the above as "C2030" format
% Another storage(such as ITSG):C20,C21,S21,C22,S22,C30,C31,S31,C32,S32,C33,S33,....
% We define the above as "C2021" format
index=1;
lm=zeros(1,2);
for i=2:lmax
    for j=0:i
        if j==0
            lm(index,:)=[i,j];
            index=index+1;
        else
            lm(index:index+1,:)=[i,j;i,j];
            index=index+2;
        end
    end
end
C2021toC2030_index=[find(lm(:,2)==0);find(lm(:,2)~=0)];
lm_=zeros(1,2);index_=1;
for i=2:lmax
    for j=0:i
        if j==0
            lm_(index_,:)=[i,j];
            index_=index_+1;
        end
    end
end
for i=2:lmax
    for j=0:i
        if j~=0
            lm_(index_:index_+1,:)=[i,j;i,j];
            index_=index_+2;
        end
    end
end
C2030toC2021_index=zeros(size(lm_,1),1);ind=1;
for ii=1:lmax
    for jj=1:size(lm_,1) 
        if lm_(jj,1)==ii  
            C2030toC2021_index(ind,1)=jj;
            ind=ind+1;
        end        
    end  
end
end

2、“Tran_vector ”函数  功能：将方阵存储的球谐系数展开为列向量

function [CS] = Tran_vector(cs,lmax,num,index)
CS=zeros((lmax+1)^2,num);
if ndims(cs)==2
    sc=gmt_cs2sc(cs);
    y_=sc2SHCvector(sc);%Do not contain C20 and degree_1 %球谐系数存放为s、c交替出现
    CS(5:(lmax+1)^2,:)=y_(index,1);%将其转换由c递增转换为s递增
    CS(1:2,:)=sc(1:2,lmax+1);CS(3,:)=sc(2,lmax+2);CS(4,:)=sc(2,lmax);
else   
for i=1:num
    sc=gmt_cs2sc(reshape(cs(i,:,:),lmax+1,lmax+1));
    y_=sc2SHCvector(sc);%Do not contain C20 and degree_1 %球谐系数存放为s、c交替出现
    CS(5:(lmax+1)^2,i)=y_(index,1);%将其转换由c递增转换为s递增
    CS(1:2,i)=sc(1:2,lmax+1);CS(3,i)=sc(2,lmax+2);CS(4,i)=sc(2,lmax);
end
end
end

核心计算和绘图程序：

clearvars -except
addpath E:\Data\result\grid
load grid_CSR.mat;load grid300_CSR.mat;load gridP4M6_CSR.mat;load grid300P4M6_CSR.mat;

name{1,1}='_CSR';name{2,1}='300_CSR';
name{3,1}='P4M6_CSR';name{4,1}='300P4M6_CSR';
k=size(name,1); num=size(grid_CSR,3);

for jj=1:k
cs=gmt_grid2cs(eval(['grid' name{jj,1};]),60);  %%转换为质量系数
for ii=1:num
    temp(:,:)=cs(ii,:,:);
    cs_temp(ii,:,:)=gmt_mc2gc(temp);   %%质量系数转换为球谐系数
end

switch jj
    case 1
        cs_CSR=cs_temp;
    case 2
        cs300_CSR=cs_temp;
    case 3
        csP4M6_CSR=cs_temp;
    otherwise
        cs300P4M6_CSR=cs_temp;
end
end

lmax=60;  [~,C2030toC2021_index]=GenerateC2030C2021TransIndex(lmax);
cs1 = Tran_vector(cs_CSR,lmax,num,C2030toC2021_index);  %%将61*61方阵的球谐展开为列向量
cs2 = Tran_vector(cs300_CSR,lmax,num,C2030toC2021_index);
cs3 = Tran_vector(csP4M6_CSR,lmax,num,C2030toC2021_index);
cs4 = Tran_vector(cs300P4M6_CSR,lmax,num,C2030toC2021_index);
csm(:,1)=mean(cs1,2);csm(:,2)=mean(cs2,2);  csm(:,3)=mean(cs3,2);csm(:,4)=mean(cs4,2);


for j=1:k
index=5;
   for i=2:lmax
     Temp(i,j)=6371e3*sqrt(mean(csm(index:index+2*i,j).^2));
     index=index+2*i+1;
   end     
end

close all;
l=2:lmax;th=1.5;
colorS=[0 1 0;0 0 1;1 0 1;0.1 0.1 0.5];

figure(1);
for i=1:k
plot(l,abs(Temp(2:end,i)),'color',colorS(i,:),'Linewidth',th);hold on;
end
grid on;
set(gca,'yscale','log');set(gca,'Fontsize',10,'ylim',[5e-8,1e-3]);
xticks([0 10 20 30 40 50 60]);
xlabel('Degree','Fontname','Time New Roman');
ylabel('Geoid degree variance','Fontname','Time New Roman');
l1=legend('Unfiltered','G300','P4M6','G300+P4M6');
set(l1,'box','off',"FontSize",10);

结果图谱：

 欢迎评论区或私信交流，多多点赞收藏，感谢支持。

参考文献：

Wahr,John.Time-variable gravity from GRACE: First results[J].Geophysical Research Letters, 2004, 31(11):293-317.DOI:10.1029/2004GL019779.