速度场模型

空间大地测量学作业学院：测绘科学与技术学院专业：测绘工程姓名：曹延虎学号：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个数据，如下图红线所示。koioinioioiiyyxxyoixoiyxQyxyxQyxv])()[(),,,(),,,(),(2221PvQPQQvQVQvTT1)(图1δ=0.01效果图图2δ=0.1效果图图3加入断层效果图1.1.1Matlab源程序（原始数据及中间值见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];fori2=1:how2forj2=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)];fork2=1:size(Unb2,1)form2=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),(),(),(jYiXhXYYXfYeXdXYcYbXaYXFfYeXdXYcYbXaYXFcYbXaYXF本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。图4拟合效果图1.2.1Matlab源程序（原始数据及中间值见EXCEL表）%40个结点局部多项式拟合的预计模型clear,clc;Y=xlsread('原始数据经纬度.xls',1,'B2:B41');%读入已知点数据,读入经度YX=xlsread('原始数据经纬度.xls',1,'C2:C41');%读入纬度X%读入40个结点数据Vy=xlsread('原始数据经纬度.xls',1,'D2:D41');%读入东西方向的速度VyVx=xlsread('原始数据经纬度.xls',1,'E2:E41');%读入南北方向的速度Vxhow3=size(Y,1);array3=zeros(how3,10);fori3=1:how3forj3=1:10ifj3==1array3(i3,j3)=1;%1elseifj3==2array3(i3,j3)=X(i3);%Xelseifj3==3array3(i3,j3)=Y(i3);%Yelseifj3==4array3(i3,j3)=X(i3)*Y(i3);%X*Yelseifj3==5array3(i3,j3)=X(i3)^2;%X^2elseifj3==6array3(i3,j3)=Y(i3)^2;%Y^2elseifj3==7array3(i3,j3)=X(i3)^2*Y(i3);%X*Y^2elseifj3==8array3(i3,j3)=Y(i3)^2*X(i3);%Y*X^2elseifj3==9array3(i3,j3)=X(i3)^3;%X^3elseifj3==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);fork=1:size(Unb3_Y)forn=1:10ifn==1Oritation(k,n)=1;%1elseifn==2Oritation(k,n)=Unl3_X(k);%Xelseifn==3Oritation(k,n)=Unb3_Y(k);%Yelseifn==4Oritation(k,n)=Unl3_X(k)*Unb3_Y(k);%X*Yelseifn==5Oritation(k,n)=Unl3_X(k)^2;%X^2elseifn==6Oritation(k,n)=Unb3_Y(k)^2;%Y^2elseifn==7Oritation(k,n)=Unl3_X(k)^2*Unb3_Y(k);%X*Y^2elseifn==8Oritation(k,n)=Unb3_Y(k)^2*Unl3_X(k);%Y*X^2elseifn==9Oritation(k,n)=Unl3_X(k)^3;%X^3elseifn==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反距离加权法拟合水平速度场反距离加权法是一个加权平均插值法，又称为谢别德法。利用此方法可以进行圆滑的或者较准确的插值。设平面上分布一系列离散点P(x,y,z)，已知其坐标为P(xi,yi)，插值zi(i=1,2,…,n)，根据周围离散点的值，通过距离加权插值求得P点的插值。nikinikiipyxdyxdzZ11)],([1)],([其中22)()(),(iiiyyxxyxd，表示由离散点到待插值点的距离。计算格网点数据时，给每个观测点赋予一个权重，所有点的权重之和为1.0。.如果某