GRACE数据的广义三角帽（TCH）计算不确定度

1、首先整理好需要计算不确定度的GRACE数据，这里以CSR level-2截断至60阶的球谐系数、DDK3和4滤波解、CSR mascon以及JPL mascon数据为例（以JPL mascon数据为参照）

y=[csr60-jpl_mas,ddk3-jpl_mas,...
   ddk4-jpl_mas,csr_mas-jpl_mas];
S=cov(y);%求取原始序列y的协方差矩阵，矩阵大小一般为n-1*n-1
n=size(S,1)+1;
save E:\RSE\result\S.mat S
J=[eye(n-1),-1*ones(n-1,1)];K=(det(S))^(1/(n-1));

2、根据TCH中的相关公式，利用迭代初值计算噪声协方差矩阵R(n*n)中的相关元素

R=zeros(n,n);
R(n,n)=1/(2* sum(sum(inv(S))));
for i=1:n-1
    for j=1:n-1
        R(i,j)=S(i,j)-R(n,n)+R(i,n)+R(j,n);
    end
end

3、matlab中根据实际情况构造目标函数

function f=fun(x)
load E:\RSE\result\S.mat;
f=(S(1,2)-x(5)+x(1)+x(2)).^2+(S(1,3)-x(5)+x(1)+x(3)).^2+...
    (S(1,4)-x(5)+x(1)+x(4)).^2+(S(2,3)-x(5)+x(2)+x(3)).^2+...
    (S(2,4)-x(5)+x(2)+x(4)).^2+(S(3,4)-x(5)+x(3)+x(4)).^2+...
    x(1).^2+x(2).^2+x(3).^2+x(4).^2;
end

4、根据Vagner G. Ferreira（2016）论文中的公式（9）构造约束条件

function [c,ceq] = mycon(x)
load E:\RSE\result\S.mat;
S_inv=inv(S);
b1=x(5)-x(1);
b2=x(5)-x(2);
b3=x(5)-x(3);
b4=x(5)-x(4);
c =b1*(S_inv(1,1)*b1+S_inv(2,1)*b2+S_inv(3,1)*b3+S_inv(4,1)*b4)+...
    b2*(S_inv(1,2)*b1+S_inv(2,2)*b2+S_inv(3,2)*b3+S_inv(4,2)*b4)+...
    b3*(S_inv(1,3)*b1+S_inv(2,3)*b2+S_inv(3,3)*b3+S_inv(4,3)*b4)+...
    b4*(S_inv(1,4)*b1+S_inv(2,4)*b2+S_inv(3,4)*b3+S_inv(4,4)*b4)-x(5);
ceq = [];
end

5、利用matlab中求解非线性多元函数最小值的fmincon函数求解在满足约束条件的情况下使得目标函数最小的参数

x0=[zeros(1,n-1),R(n,n)]';
A=[];   
Aeq=[];
b=[];
beq=[];
lb=zeros(1,n);
ub=[];
[x,fval] = fmincon('fun',x0,[],[],[],[],lb,[],'mycon')
%fun为目标函数，mycon为约束条件，求解参数x的下限为0
R=zeros(n,n);
R(1,n)=x(1,1);R(2,n)=x(2,1);R(3,n)=x(3,1);R(4,n)=x(4,1);R(5,n)=x(5,1);
R(n,1)=x(1,1);R(n,2)=x(2,1);R(n,3)=x(3,1);R(n,4)=x(4,1);
for i=1:n-1
    for j=1:n-1
        R(i,j)=S(i,j)-R(n,n)+R(i,n)+R(j,n);
    end
end
det(R)
%为验证结果的正确性，观察矩阵R的行列式是否大于0

 参考文献

[1] Tavella P ,  Premoli A . Estimating the Instabilities of N Clocks by Measuring Differences of their Readings[J]. Metrologia, 1994, 30(5):479.
[2]姚朝龙, 李琼, 罗志才,等. 利用广义三角帽方法评估GRACE反演中国大陆地区水储量变化的不确定性[J]. 地球物理学报, 2019, 62(3):15.
[3] Ferreira V G ,  Montecino H ,  Yakubu C I , et al. Uncertainties of the Gravity Recovery and Climate Experiment time-variable gravity-field solutions based on three-cornered hat method[J]. Journal of Applied Remote Sensing, 2016, 10(1):015015.