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.