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-1for j=1:n-1R(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-1for j=1:n-1R(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.

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

pandas使用groupby函数、agg函数获取每个分组聚合对应的均值（mean）实战：计算分组聚合单数据列的均值、计算分组聚合多数据列的均值
pandas使用groupby函数.agg函数获取每个分组聚合对应的均值(mean)实战:计算分组聚合单数据列的均值.计算分组聚合多数据列的均值目录
SQL Server中的报表–如何使用数据透视表和日期计算来获取有价值的报表
介绍 (Introduction) A few months back I had been working on an interesting proof of concept for a huma ...
云原生数据湖以存储、计算、数据管理等能力通过信通院评测认证
又一项大能力-云原生数据湖获得信通院认证啦! 近日,中国信息通信研究院 (以下简称"信通院") 正式公布了第十四批"大数据产品能力评测"结果,腾讯云云原生数据湖 ...
三种数据交换方式的时延计算
三种数据交换方式的时延计算 part1:什么是时延?有哪三种数据交换方式? 时延指的是计算机网络性能的一种,表示数据从网络的一端传送到另一端所用的时间. 三种数据交换方式:电路交换.报文交换.分组交换 ...
使用Python对股票数据进行数据分析(一）-计算日线行情、5日均线、10日均线行情并显示
使用Python对股票数据进行数据分析(一)-计算日线行情.5日均线.10日均线行情并显示各种炒股软件上可以显示各种技术指标,可以帮助投资者进行技术分析.这些股市中的这些指标都是怎么计算出来的呢?这 ...
计算机表格中如何计算数据透视表,如何在EXCEL数据透视表中进行计算 |
excel 数据透视表中如何插入公式数据透视表>公式>计算字段 Excel 数据透视表中如何算占比? 在表格中右键透视表之后,打开值字置,点击值显示方式,然下拉菜单择多种占比. ex ...
计算机表格中如何计算数据透视表,Excel表格中在数据透视表中添加计算字段的方法...
计算字段是使用数据透视表中的字段同其他内容经过计算后得到的,如果用户需要在数据透视表中自定义计算公式以计算数据,可以通过添加计算字段来实现,下面介绍Excel表格中在数据透视表中添加计算字段的具体操作 ...
计算机表格中如何计算数据透视表,Excel中如何在数据透视表中进行计算
会计工作中离不开excel电子表格软件,它不仅具有数据输入.输出.显示.分类.统计.查询等数据处理的基本功能,还具有强大的数据分析功能与程序执行自动化功能,为会计人员的工作提供了许多便利.数据透视表是 ...
【个推CTO谈数据智能】之数据安全计算体系
作者|个推CTO 安森引言本文是数据智能系列的第四篇.前三篇文章(<数据智能时代来临:本质及技术体系要求><多维度分析系统的选型方法> <我们理解的数据中台> ...

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

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

最新文章

热门文章