利用GRACE球谐数据计算地表位移的基本原理与实现

1.基本理论

由于地表质量发生变化导致地表的负荷变形，利用GRACE数据计算地表的垂直（2）、水平东向（3）和水平北向位移（4）方法如下：

我们之前已经利用冯伟老师提供的MATLAB工具包计算了EWH，下面我们具体看看两类计算公式的异同点：

我们可以发现，有一些部分是一致的，即

唯一不同的就是前面乘以的系数。我们看的，垂直位移的系数是

a*h/(1+k)

而等效水高形式的系数是

ap(2l+1)/3pw(k+1)

因此只要在冯伟老师的代码中修改这部分系数就可以利用GRACE球谐系数计算垂直位移。

然而对于水平位移涉及到勒让德函数的偏导数，下面的公式给出了计算勒让德偏一阶和二阶偏导数的递推公式：

2.计算勒让德函数的程序

我们先给出常用的勒让德函数的递推公式，然后再给出勒让德函数求偏导数的程序

勒让德函数递推公式：

function le=legendre_n(x,Nmax)
% 2024-3-20
le(1:Nmax+1,1:Nmax+1)=0;

% l = m =0
le(1,1)=1.0;
% l = 1,m = 0         
le(2,1)=sqrt(3.0)*x;
% l = m =1         
le(2,2)=sqrt(3.0)*sqrt(1.0-x*x);
% l = m
for i=3:Nmax+1
    le(i,i)=sqrt( (1.0-x*x) * (2*i-1) / (2*i-2) )*le(i-1,i-1);
end
% n-m = 1
for i=3:Nmax+1
    le(i,i-1)=x*sqrt(2*i-1)*le(i-1,i-1);
end
%others
for j=1:Nmax-1
    for i=j+2:Nmax+1
        le(i,j)=x*sqrt((2.0*i-3)*(2.0*i-1)/((i-j)*(i+j-2)))*le(i-1,j) ...
            -sqrt((2.0*i-1)*(i+j-3.0)*(i-j-1.0)/((2.0*i-5)*(i+j-2)*(i-j)))*le(i-2,j);
    end
end
 

勒让德函数一阶偏导：

function le_dtheta=legendre_partial_n(x,Nmax)
% 2024-3-20
rP = legendre_n(x,Nmax);

le_dtheta = zeros(Nmax+1);
for i = 0:Nmax
    for j = 0:i
        if i == j
            le_dtheta(i+1,i+1) = i*x/sqrt(1-x^2) * rP(i+1,i+1);
        else
            le_dtheta(i+1,j+1) = i*x/sqrt(1-x^2) * rP(i+1,j+1)...
                -1/sqrt(1-x^2)*sqrt( (i-j)*(i+j)*(2*i+1)/(2*i-1) ) * rP(i,j+1);
        end
    end
end
                
 

勒让德函数二阶偏导：

function le_dtheta=legendre_partial_nn(x,Nmax)
% 2024-3-20

rP = legendre_n(x,Nmax);

le_dtheta = zeros(Nmax+1);
for i = 0:Nmax
    for j = 0:i
        if i == j
            le_dtheta(i+1,i+1) = i*i-j-j*j*x/(1-x^2) * rP(i+1,i+1);
        else
            le_dtheta(i+1,j+1) = i*i-j-j*j*x/(1-x^2) * rP(i+1,j+1)...
                +x/(1-x^2)*sqrt( (i-j)*(i+j)*(2*i+1)/(2*i-1) ) * rP(i,j+1);
        end
    end
end
                
 

实际运算例子：

lon = 0.5:1:359.5;lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);

Nmax = 60; % max degree
Nlat = size(lon,1); % size of lat grid

for ilat = 1: Nlat
    rCoLat= 90.0-lat(ilat,1);
    rCoLat2= cos(deg2rad(rCoLat));
    rP = legendre_n( rCoLat2, Nmax );
    rP_d1 = legendre_partial_n( rCoLat2, Nmax );
    rP_d2 = legendre_partial_nn( rCoLat2, Nmax );
end

可以得到61×61阶的结果，这与60阶GRACE球谐数据匹配。

3.计算GRACE数据得到的地表负荷位移

这次之前，我们需要弄清楚一个事实，即大地水准面系数和质量系数之间的关系。

回到Whar et al.(1998)年的论文，下式表达了两者之间的关系，即我们处理的GRACE数据是大地水准面系数，而利用GLDAS等水文模型导出的是质量球谐系数。

下面我们利用GLDAS数据计算地表位移。利用冯伟老师的工具箱，我们可以得到球谐展开的GLDAS系数，该系数是质量系数，即上式的右边，而如果要计算地表位移，需要首先将其转成大地水准面系数（左边），然后计算负荷位移。上一篇专栏我介绍了一个修改的，可以将GLDAS水文模型球形展开的例子，感兴趣的可以阅读一下。在冯伟老师的工具箱中，具体需要修改的函数是gmt_gc2mc。这里我编写一个gmt_mc2gc的函数，结合gmt_gc2lc即可得到负荷垂直位移。

function outcoeff = gmt_mc2gc(incoeff)

% Converts geoid coefficients (gc) to mass coefficients (mc)
% FENG Wei 19/04/2015
% State Key Laboratory of Geodesy and Earth's Dynamics
% Institute of Geodesy and Geophysics, Chinese Academy of Sciences
% [email protected]
% modified 20/3/2024 Chistrong Wen
[rows,cols] = size(incoeff);
if rows == cols                    % field is in CS-format
    maxdeg  = rows - 1;
    %    incoeff = cs2sc(incoeff);
elseif cols-2*rows == -1            % field is in SC-format already
    maxdeg  = rows - 1;
    incoeff = gmt_sc2cs(incoeff);            % convert to CS-format
else
    error('Check format of field.')
end

ae        =  6378136.3;             % semi-major openaxis of ellipsoid [m]
rho_w = 1000;       % density of water kg/m^3
rho_ave = 5517;     % average density of the earth kg/m^3

%% mc2load
dens = (3.*rho_w)/(ae*rho_ave);
% Han and Wahr Love numbers
load LoveNumbers_lln.mat
% %Refer to Chamnbers deg1_coef.txt  Fengwei
k(2) = 0.021;
sc = gmt_cs2sc(incoeff);
% loop over degrees to multiply each row of coefficients in sc-format
% with the factor from equation (13) of Wahr et al. 1998
for ll = 0:1:maxdeg
    factor(ll+1)  = dens*( 1+k(ll+1) )/( 2*ll+1 );
    scnew(ll+1,:) = factor(ll+1)*sc(ll+1,:);
end

outcoeff = gmt_sc2cs(scnew);
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CCBY Chistrong Wen
% University of Chinese Academy of Sciences
% 2024-3-20
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GLDAS INFO
O.lon = mascon_csr.lon;O.lat = mascon_csr.lat;O.rg = mascon_csr.rg(:,:,12);
figure
subplot(3,1,1),wzq_plot(O),title('surface mass'),caxis([-20,20]),colorbar
%% expand to 60 degree
grid = O.rg;
degree_max = 60;
[cs] = gmt_grid2cs(grid,degree_max);
%% applying filter method
destrip_method = 'NONE';
radius_filter = 000; 
type = 1;

gc_cs = gmt_mc2gc(cs);
load_cs = gmt_mc2load(gc_cs);

outcoeff = gmt_gc2lc(gc_cs);

[grid_filter_load]=gmt_cs2grid(outcoeff,radius_filter,type,destrip_method);
[grid_filter_mass]=gmt_cs2grid(cs,radius_filter,type,destrip_method);

lon = 0.5:1:359.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
subplot(3,1,3)
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg  = grid_filter_load;
wzq_plot(out_soil_filter)
title('degree 60 mass'),colorbar

subplot(3,1,2)
figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg  = grid_filter_mass;
wzq_plot(out_soil_filter)
title('mass load'),colorbar

结果图：

代码包见资源。

♥欢迎点赞收藏♥