引言
由于GPS观测点分布离散且不均匀,在进行应变计算和分析前,需要对速度场进行插值,获得均匀分布的速度场。一般采用Kriging 法估计在均匀网格节点上的速度值,需要下载克里金MATLAB工具箱(matlab_Kriging_toolbox)并添加到路径中 (下载链接:https://pan.baidu.com/s/1p0lt2G8KQ-els7Dyw5jPFg,提取码:wcss)
克里金插值原理
网上介绍克里金插值的帖子和文章非常多,这里推荐几个供大家参考:
-
对理解公式原理非常有帮助的:
https://xg1990.com/blog/archives/222
https://zhuanlan.zhihu.com/p/89619040 -
比较详细介绍Dace中例子的文章:https://blog.csdn.net/qq_40937675/article/details/89792122
克里金插值步骤
- (1)对于观测数据,两两计算距离与半方差
- (2)寻找一个拟合曲线拟合距离与半方差的关系,从而能根据任意距离计算出相应的半方差
- (3)计算出所有已知点之间的半方差 rij
- (4)对于未知点 zo,计算它到所有已知点 zi的半方差 rio
- (5)求解第四节中的方程组,得到最优系数 λi 使用最优系数
- (6)对已知点的属性值进行加权求和,得到未知点 zo的估计值
实例阐述
在具体用matlab实现时分为以下几步:
- 1)导入观测数据文件,包含点位坐标和观测值,建立变量。
在GNSS观测值中包含北东高三个方向的形变量,由于GNSS观测对于垂直向形变不敏感,因此在插值的时候是使用了水平两个方向的形变量。 - 2)设置范围和格网大小使用工具箱进行插值。
- 3)根据点位误差大小调整拟合函数及其他参数。
Dace工具箱中对于回归函数包含0阶、1阶和2阶多项式三种,对于相关函数包含指数函数、高斯函数、线性函数、球谐函数等六种,每种函数的拟合效果不同,可以根据实际的数据拟合情况来酌情选择合适的函数。
- 4)输出插值点文件并进行可视化,为便于绘制各点的误差圆,我是输出文件后在GMT中进行可视化。
附件代码如下,可先尝试Dace工具箱中的例子,之后按照其例子修改:
%% 导入文件,切出各列
GPS_file=xlsread('GPSdata.xlsx');
lon=GPS_file(:,1);
lat=GPS_file(:,2);
v_e=GPS_file(:,3);
v_n=GPS_file(:,4);
v_h=sqrt(v_e.*v_e+v_n.*v_n);
sigma_e=GPS_file(:,<