速度场模型

空间大地测量学作业

学院：测绘科学与技术学院

专业：测绘工程

姓名：曹延虎

学号：201410554

2015年7月

地壳运动速度场模型的建立方法

一 模型实验分析

1.1 多面函数法拟合水平速度场

多面函数法是美国Hardy 教授于1977年提出的，它基于任何一个光滑的数学表面总可用一系列规则的数学函数以任意精度逼近的思想，利用已测点推估未测点。由于具有设计灵活、可控性强等优点，多面函数法自提出以来，就被广泛用于与地学有关的插值问题。

设有m 个已测点s(x,y)，s 为非随机信号，s(x,y)为点(x,y)上的水平运动速率，可用n 个核函数的总和去逼近任意点的水平运动速率。式中ｎ为所选的已知水平运动速率的点，称为结点；αi 为待定参数；Q 为核函数，一般选为对称距离型函数，δ为平滑因子。

本次模型取k=1/2，n=20，δ=0.01和δ=0.1的正双曲面函数为核函数。平滑因子δ的作用是改变核函数的形状。核函数一旦选定，平滑因子的改变同样也会影响内插与拟合的效果。一般来讲，δ越大，核函数所表达的曲面越平缓；δ越小，曲面越陡峭。不同的核函数对平滑因子的敏感度不同，显然，平滑因子的敏感度越低越好，这样的拟合效果将更趋于稳定。 本次实验拟合了30个数据，如下图红线所示。 k oi oi n

i oi oi i y y x x yoi xoi y x Q y x y x Q y x v ])()[(),,,()

,,,(),(2221δα+-+-=⋅=∑

=Pv Q PQ Q v Q V Q v T T 1)(-=-=⋅=ααα

图1 δ=0.01 效果图

图2 δ=0.1 效果图

图3 加入断层效果图

1.1.1 Matlab源程序（原始数据及中间值见EXCEL表）

%读入40个已知数据

Jb2=xlsread('原始数据经纬度.xls',1,'B2:B41'); %读入经度

Jl2=xlsread('原始数据经纬度.xls',1,'C2:C41'); %读入纬度

%读入20个结点数据

JJb2=xlsread('原始数据经纬度.xls',1,'I2:I21');%读入经度

JJl2=xlsread('原始数据经纬度.xls',1,'J2:J21');%读入纬度

how2=size(Jb2,1);

vol2=size(JJb2,1);

array2=[how2,vol2];

for i2=1:how2

for j2=1:vol2

array2(i2,j2)=sqrt((Jb2(i2)-JJb2(j2))^2+(Jl2(i2)-JJl2(j2))^2+0.01^2); end

end

%读入已知点的速度值

Vx2=xlsread('原始数据经纬度.xls',1,'D2:D41');

Vy2=xlsread('原始数据经纬度.xls',1,'E2:E41');

%计算拟合模型的参数

Rx2=inv(array2'*array2)*array2'*Vx2;

Ry2=inv(array2'*array2)*array2'*Vy2;

xlswrite('原始数据经纬度.xls',Rx2,2,'A1');%将X方向的参数值输出

xlswrite('原始数据经纬度.xls',Ry2,2,'B1');%将Y方向的参数值输出

%计算改正数VVx2和VVy2

VVx2=array2*Rx2-Vx2;

VVy2=array2*Ry2-Vy2;

%计算单位权中误差

dx2=sqrt((VVx2'*VVx2)/(how2-vol2));

dy2=sqrt((VVy2'*VVy2)/(how2-vol2));

xlswrite('原始数据经纬度.xls',dx2,2,'C1'); %将X方向的中误差输入到excle中去xlswrite('原始数据经纬度.xls',dy2,2,'D1'); %将Y方向的中误差输入到excle中去

%计算待估参数的协因数,两个协因数阵相同

Qx2=inv(array2'*array2);

Qy2=inv(array2'*array2);

Dx2=sqrt(Qx2)*dx2;%计算方差

Dy2=sqrt(Qy2)*dy2;

%将参数的中误差输出

xlswrite('原始数据经纬度.xls',Dx2,3,'E1'); %将X方向的参数中误差输入到excle中去xlswrite('原始数据经纬度.xls',Dy2,3,'E41'); %将Y方向的参数中误差输入到excle中去

%用模型计算待定点的速度

Unb2=xlsread('原始数据经纬度.xls',1,'D127:D157');%读入待求点的经度

Unl2=xlsread('原始数据经纬度.xls',1,'E127:E157');%读入待求点的纬度

Array_Un2=[size(Unb2,1),size(JJb2,1)];

for k2=1:size(Unb2,1)

for m2=1:size(JJb2,1)

Array_Un2(k2,m2)=sqrt((Jb2(k2)-JJb2(m2))^2+(Jl2(k2)-JJl2(m2))^2+0.01^2);

end

end

Speed_X2=Array_Un2*Rx2; %计算南北方向（X方向）的速度

Speed_Y2=Array_Un2*Ry2; %计算东西方向（Y方向）的速度

%往excle中写数据

%xlswrite(filename,A,sheet,range); % A就是待写的数据

xlswrite('原始数据经纬度.xls',Speed_X2,1,'H127');

xlswrite('原始数据经纬度.xls',Speed_Y2,1,'I127');

1.2 局部多项式法拟合水平速度场

多项式的形式有三种：一次多项式、二次多项式和三次多项式。在插值时，寻找一个适合的函数并非易事，如多项式阶数太大的问题，其波动也会很大。针对这个问题，局部多项式法是一个不错的选择。它作为常用的拟合方法之一，能在给定搜索区域内对插值对象所有点插值出合理特定阶数的多项式，局部多项式插值产生的曲面根多依赖于局部的变异。常见的多项式形式如下：