一、數(shù)據(jù)準(zhǔn)備

首先將從ICGEM官網(wǎng)（http://icgem.gfz-potsdam.de/home）下載的GRACE(-FO)數(shù)據(jù)解壓存至一個(gè)文件夾下：

二、數(shù)據(jù)格式介紹

從ICGEM下載的GRACE（-FO）文件有著相同的數(shù)據(jù)格式，其文件名如“GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gfc”, 其中“GSM”表示該文件對(duì)應(yīng)的是大地水準(zhǔn)面球諧位系數(shù)，“02095-2002120”該數(shù)據(jù)文件對(duì)應(yīng)的數(shù)據(jù)觀測(cè)時(shí)間是2002年第95天至2002年第120天，“UTCSR”表示的是該文件由CSR機(jī)構(gòu)發(fā)布。

打開(kāi)任意一期數(shù)據(jù)文件，可發(fā)現(xiàn)數(shù)據(jù)由頭文件、球諧系數(shù)組成:

頭文件介紹了GRACE數(shù)據(jù)的技術(shù)說(shuō)明文檔以及其他的基本信息，此處不再進(jìn)行詳細(xì)闡述。紅色框中闡述了一些基本常數(shù)，如地球半徑、地球密度常數(shù)、球諧系數(shù)文件的最大階數(shù)等。最后是極為重要的Stocks球諧系數(shù)，主要包含7列數(shù)據(jù)：第一列無(wú)具體含義、第二列L對(duì)應(yīng)階數(shù)、第三列M對(duì)應(yīng)次數(shù)、第四列C對(duì)應(yīng)Cnm球諧系數(shù)、第五列S對(duì)應(yīng)Snm球諧系數(shù)，第六列sigma C對(duì)應(yīng)Cnm的誤差，第7列sigma S對(duì)應(yīng)Snm的誤差。

指定讀取球諧系數(shù)的最大階數(shù)，如maxDegree=60,并指定讀取的數(shù)據(jù)時(shí)間范圍，如timegap=[200204 202012]，將下載的GRACE（-FO）數(shù)據(jù)文件夾的路徑賦給dir_in，如“dir_in='H:\matlab_text_code\GRACE_and_GRACE_Fo\CSR_06_Grace_and_Grace_Fo'；”，運(yùn)行[GCSR_cnm,GCSR_snm,GCSR_time,GCSR_time_ym]=hfl_readgsm(dir_in,'ICGEM',maxDegree,timegap);即可獲取對(duì)應(yīng)時(shí)間段的GRACE（-FO）球諧系數(shù)，以及對(duì)應(yīng)的時(shí)間。

附上hfl_readgsm.m和hfl_readstocks.m函數(shù)：

function [cnm,snm,time,time_ym]=hfl_readgsm(dir_in,type,lmax,timegap)

?

% Read the gravity filed files

%

% INPUT:

%?? dir_in????? full path

%

% OUTPUT:

%? cnm&snm? spherical harmonic coefficients|c\| and |s\|

%?? time??????? fraction in the year

%?? time_ym?? year+month

%

?

switch (type)

??? case 'ICGEM'

??????? files=dir(fullfile(dir_in,'*.gfc'));

??????? lenfiles=length(files);

??????? cnm=zeros(lmax+1,lmax+1,lenfiles);snm=cnm;time_bounds=zeros(lenfiles,2);

??????? h = waitbar(0,'Loading the GRCAE data...');

??????? for i=1:lenfiles

??????????? filename=fullfile(dir_in,files(i).name);

??????????? time_bounds(i,:)=[str2double(files(i).name(7:13)) str2double(files(i).name(15:21))];

??????????? [cnm(:,:,i),snm(:,:,i)]=hfl_readstocks(filename,lmax);

??????????? waitbar(i/lenfiles,h,['loading' num2str(i) '/ ' num2str(lenfiles)]);

??????? end

??????? close(h)

??????? [time,time_ym]=gmt_TimeBoundsToGRACETime(time_bounds);

??????? pos=(time_ym>=timegap(1)&time_ym<=timegap(2));

??????? cnm=cnm(:,:,pos);

??????? snm=snm(:,:,pos);

??????? time=time(pos);

??????? time_ym=time_ym(pos);

??? otherwise

??????? warndlg('Format of filetype is only ICGEM!','Warning');

end

end

??

?????

function [cnm,snm]=hfl_readstocks(filename,lmax)

?

% Read stocks potential coeffcients from file.

%

% Input:

% filename: file of potential coeffcients in ICGEM format (*.gfc).

% Output:

% cnm,snm: matricies containing the potential coefficients.

% The dimensions are (n + 1) x (n + 1) with n = lmax.

% The lower triangular matrix elements cnm(n + 1, m + 1) and snm(n + 1, m + 1)

% contain the potential coefficients of degree n and order m.

?

% open the file

fid = fopen(filename);

hasErrors = true;

% read header

while 1

? line = fgetl(fid);

? if ~ischar(line), break, end

? if isempty(line), continue, end

? keyword = textscan(line,'%s',1);

?

? if(strcmp(keyword{1}, 'max_degree'))

??? cells = textscan(line,'%s%d',1);

??? maxDegree = cells{2};

? end

?

? if(strcmp(keyword{1}, 'errors'))

??? cells = textscan(line,'%s%s',1);

??? if(strcmp(cells{2}, 'no'))

????? hasErrors = false;

??? end

? end

?

? if(strcmp(keyword{1}, 'end_of_head'))

??? break

? end

end

?

if (lmax==maxDegree)

???

??? % init the output matricies.

??? cnm = zeros(lmax+1, lmax+1);

??? snm = zeros(lmax+1, lmax+1);

???

??? % read potential coefficients

??? while 1

??????? line = fgetl(fid);

??????? if ~ischar(line), break, end

??????? if isempty(line), continue, end

??????? if(hasErrors)

??????????? cells = textscan(line,'%s%d%d%f%f%f%f');

??????? else

??????????? cells = textscan(line,'%s%d%d%f%f');

??????? end

??????? if(cells{2}+1<=lmax+1&&cells{3}+1<=lmax+1)

??????????? cnm(cells{2}+1, cells{3}+1) = cells{4};

??????????? snm(cells{2}+1, cells{3}+1) = cells{5};

??????? else

?????? ?????continue;

??????? end

??? end

else

??? warndlg('Format of file is wrong, please check the maxDegree!','Warning');

end

fclose(fid);

end

???