速度场模型

空间大地测量学作业学院：测绘科学与技术学院专业：测绘工程姓名：曹延虎学号：2014105542015年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 ni 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:how2for j2=1:vol2array2(i2,j2)=sqrt((Jb2(i2)-JJb2(j2))^2+(Jl2(i2)-JJl2(j2))^2+0.01^2); endend%读入已知点的速度值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和VVy2VVx2=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);endendSpeed_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 局部多项式法拟合水平速度场多项式的形式有三种：一次多项式、二次多项式和三次多项式。

在插值时，寻找一个适合的函数并非易事，如多项式阶数太大的问题，其波动也会很大。

针对这个问题，局部多项式法是一个不错的选择。

它作为常用的拟合方法之一，能在给定搜索区域内对插值对象所有点插值出合理特定阶数的多项式，局部多项式插值产生的曲面根多依赖于局部的变异。

常见的多项式形式如下：33222222),(),(),(jY iX hXY Y X fY eX dXY cY bX a Y X F fY eX dXY cY bX a Y X F cYbX a Y X F +++++++++=+++++=++=本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。

图4 拟合效果图1.2.1 Matlab 源程序（原始数据及中间值见EXCEL 表）%40个结点局部多项式拟合的预计模型clear,clc;Y=xlsread('原始数据经纬度.xls',1,'B2:B41'); %读入已知点数据,读入经度 Y X=xlsread('原始数据经纬度.xls',1,'C2:C41'); %读入纬度 X %读入40个结点数据Vy=xlsread('原始数据经纬度.xls',1,'D2:D41'); %读入东西方向的速度 Vy Vx=xlsread('原始数据经纬度.xls',1,'E2:E41'); %读入南北方向的速度 Vx how3=size(Y,1);array3=zeros(how3,10);for i3=1:how3for j3=1:10if j3==1array3(i3,j3)=1; % 1elseif j3==2array3(i3,j3)=X(i3); % Xelseif j3==3array3(i3,j3)=Y(i3); % Yelseif j3==4array3(i3,j3)=X(i3)*Y(i3); % X*Yelseif j3==5array3(i3,j3)=X(i3)^2; %X^2elseif j3==6array3(i3,j3)=Y(i3)^2; %Y^2elseif j3==7array3(i3,j3)=X(i3)^2*Y(i3); % X*Y^2elseif j3==8array3(i3,j3)=Y(i3)^2*X(i3); % Y*X^2elseif j3==9array3(i3,j3)=X(i3)^3; % X^3elseif j3==10array3(i3,j3)=Y(i3)^3; % Y^3endendend%计算拟合模型的参数Rx3=inv(array3'*array3)*array3'*Vx; %计算X方向的参数Ry3=inv(array3'*array3)*array3'*Vy; %计算Y方向的参数xlswrite('原始数据经纬度.xls',Rx3,2,'A1');%将X方向的参数值输出xlswrite('原始数据经纬度.xls',Ry3,2,'B1');%将Y方向的参数值输出%计算改正数VVx2和VVy2VVx3=array3*Rx3-Vx;VVy3=array3*Ry3-Vy;% %计算单位权中误差dx3=sqrt((VVx3'*VVx3)/(how3-10));dy3=sqrt((VVy3'*VVy3)/(how3-10));xlswrite('原始数据经纬度.xls',dx3,2,'C1'); %将X方向的中误差输入到excle中去 xlswrite('原始数据经纬度.xls',dy3,2,'D1'); %将Y方向的中误差输入到excle中去% %计算待估参数的协因数,两个协因数阵相同Qx3=inv(array3'*array3);Qy3=inv(array3'*array3);% %计算参数的中误差Dx3=sqrt(Qx3)*dx3;Dy3=sqrt(Qy3)*dy3;% %将参数的中误差输出xlswrite('原始数据经纬度.xls',Dx3,3,'A1'); %将X方向的参数中误差输入到excle 中去xlswrite('原始数据经纬度.xls',Dy3,3,'B1'); %将X方向的参数中误差输入到excle 中去% %用模型计算待定点的速度% %计算南北方向（X方向）的速度Unb3_Y=xlsread('原始数据经纬度.xls',1,'D127:D157');%读入待求点的经度 YUnl3_X=xlsread('原始数据经纬度.xls',1,'E127:E157');%读入待求点的纬度 XOritation=zeros(size(Unb3_Y),10);for k=1:size(Unb3_Y)for n=1:10if n==1Oritation(k,n)=1; % 1elseif n==2Oritation(k,n)=Unl3_X(k); % Xelseif n==3Oritation(k,n)=Unb3_Y(k); % Yelseif n==4Oritation(k,n)=Unl3_X(k)*Unb3_Y(k); % X*Yelseif n==5Oritation(k,n)=Unl3_X(k)^2; %X^2elseif n==6Oritation(k,n)=Unb3_Y(k)^2; %Y^2elseif n==7Oritation(k,n)=Unl3_X(k)^2*Unb3_Y(k); % X*Y^2elseif n==8Oritation(k,n)=Unb3_Y(k)^2*Unl3_X(k); % Y*X^2elseif n==9Oritation(k,n)=Unl3_X(k)^3; % X^3elseif n==10Oritation(k,n)=Unb3_Y(k)^3; % Y^3endendendX_Speed=Oritation*Rx3; %计算相应点的X 方向速度Y_Speed=Oritation*Ry3; %计算相应点的Y 方向速度xlswrite('原始数据经纬度.xls',X_Speed,1,'F127'); %将X 方向的参数中误差输入到excle 中去xlswrite('原始数据经纬度.xls',Y_Speed,1,'G127'); %将Y 方向的参数中误差输入到excle 中去1.3 反距离加权法拟合水平速度场反距离加权法是一个加权平均插值法，又称为谢别德法。