速度场模型.docx

1、速度场模型空 间 大 地 测 量 学 作 业 学院：测绘科学与技术学院专业： 测绘工程 姓名： 曹延虎 学号： 201410554 2015年7月地壳运动速度场模型的建立方法一 模型实验分析1.1 多面函数法拟合水平速度场多面函数法是美国Hardy教授于1977年提出的，它基于任何一个光滑的数学表面总可用一系列规则的数学函数以任意精度逼近的思想，利用已测点推估未测点。由于具有设计灵活、可控性强等优点，多面函数法自提出以来，就被广泛用于与地学有关的插值问题。 设有m个已测点s(x,y)，s为非随机信号，s(x,y)为点(x,y)上的水平运动速率，可用n个核函数的总和去逼近任意点的水平运动速率。式

2、中为所选的已知水平运动速率的点，称为结点；i为待定参数；Q为核函数，一般选为对称距离型函数，为平滑因子。本次模型取k=1/2，n=20，=0.01和=0.1的正双曲面函数为核函数。平滑因子的作用是改变核函数的形状。核函数一旦选定，平滑因子的改变同样也会影响内插与拟合的效果。一般来讲，越大，核函数所表达的曲面越平缓；越小，曲面越陡峭。不同的核函数对平滑因子的敏感度不同，显然，平滑因子的敏感度越低越好，这样的拟合效果将更趋于稳定。本次实验拟合了30个数据，如下图红线所示。 图1 =0.01 效果图图2 =0.1 效果图 图3 加入断层效果图1.1.1 Matlab源程序（原始数据及中间值见EXCE

3、L表）%读入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(

4、j2)2+(Jl2(i2)-JJl2(j2)2+0.012); 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

5、 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

6、)*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

7、=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.012); end endSpeed_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(原始数据经纬

8、度.xls,Speed_Y2,1,I127);1.2 局部多项式法拟合水平速度场多项式的形式有三种：一次多项式、二次多项式和三次多项式。在插值时，寻找一个适合的函数并非易事，如多项式阶数太大的问题，其波动也会很大。针对这个问题，局部多项式法是一个不错的选择。它作为常用的拟合方法之一，能在给定搜索区域内对插值对象所有点插值出合理特定阶数的多项式，局部多项式插值产生的曲面根多依赖于局部的变异。常见的多项式形式如下： 本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。图4 拟合效果图1.2.1 Matlab源程序（原始数据及中间值见EXCEL表）%

9、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:how3 for j3=1:10 if j3=1 a

10、rray3(i3,j3)=1; % 1 elseif j3=2 array3(i3,j3)=X(i3); % X elseif j3=3 array3(i3,j3)=Y(i3); % Y elseif j3=4 array3(i3,j3)=X(i3)*Y(i3); % X*Y elseif j3=5 array3(i3,j3)=X(i3)2; %X2 elseif j3=6 array3(i3,j3)=Y(i3)2; %Y2 elseif j3=7 array3(i3,j3)=X(i3)2*Y(i3); % X*Y2 elseif j3=8 array3(i3,j3)=Y(i3)2*X(i3)

11、; % Y*X2 elseif j3=9 array3(i3,j3)=X(i3)3; % X3 elseif j3=10 array3(i3,j3)=Y(i3)3; % Y3 end end end %计算拟合模型的参数 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和VVy2

12、 VVx3=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);% %计算参数的中误差 D

13、x3=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);%读入待求点的经度 Y Unl3_X=xlsread(原始数据经纬度.xls,1,E127:E157);%读入待求点的纬度 X Oritatio

14、n=zeros(size(Unb3_Y),10); for k=1:size(Unb3_Y) for n=1:10 if n=1 Oritation(k,n)=1; % 1 elseif n=2 Oritation(k,n)=Unl3_X(k); % X elseif n=3 Oritation(k,n)=Unb3_Y(k); % Y elseif n=4 Oritation(k,n)=Unl3_X(k)*Unb3_Y(k); % X*Y elseif n=5 Oritation(k,n)=Unl3_X(k)2; %X2 elseif n=6 Oritation(k,n)=Unb3_Y(k)2

15、; %Y2 elseif n=7 Oritation(k,n)=Unl3_X(k)2*Unb3_Y(k); % X*Y2 elseif n=8 Oritation(k,n)=Unb3_Y(k)2*Unl3_X(k); % Y*X2 elseif n=9 Oritation(k,n)=Unl3_X(k)3; % X3 elseif n=10 Oritation(k,n)=Unb3_Y(k)3; % Y3 end end end X_Speed=Oritation*Rx3; %计算相应点的X方向速度 Y_Speed=Oritation*Ry3; %计算相应点的Y方向速度 xlswrite(原始数据

16、经纬度.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点的插值。 其中，表示由离散点到待插值点的距离。计算格网点数据时，给每个观测点赋予一个权重，所有点的权重之和为1.0。

17、.如果某个观测点与格网点重合，那么这个观测点的权重为1.0，其它观测点的权重均为0.0。 该方法的一个特点就是在格网区域内，观测点周围将为产生“牛眼”现象。我们可以通过一定的方法进行消除，比如在进行格网化的时候，可以选择一个圆滑参数，该参数只要大于零，则能保证对于一个特定的节点，所有的观测点都不能被赋予权重1.0，包括该点与结点重合。本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。图5 拟合效果图1.3.1 C+源程序（原始数据见文件夹）#include #include #include #include using namespace

18、std;#define infile1 C:UserspcDesktop反距离加权法原始数据.txt /起算数据文件#define outfile1 C:UserspcDesktop反距离加权法拟合数据.txt /未知点速度场文件const int QS=40; /读取起算数据的个数const int Fix=5; /需要参考点个数int main() ifstream inf1; ofstream ouf1; inf1.open(infile1); if (inf1.fail() cout can not open infile1:endl; ouf1.open(outfile1); if

19、(ouf1.fail() cout can not open outfile1: l0 b0 Espeed0 Nspeed0; while(!inf1.eof() +i; inf1 li bi Espeedi Nspeedi; /定义未知点 double lo=0,bo=0;/计算待插值点，循环开始！for (m=0;m13;m+) for (n=0;n21;n+) lo=100+0.5*n,bo=29+0.5*m; /定义存储待插值点到起算点的距离数组 double DQS=0.0,D_BFQS=0.0; /计算待插值点到起算点的距离 for(i=0;iQS;i+) Di=1/(lo-li)

20、*(lo-li)+(bo-bi)*(bo-bi); D_BFi=Di; /对数组D进行从大到小排序 double nextElement=0; double arrayQS+1=0.0; for (i=0;iQS;i+) arrayi+1=Di; i=0,j=0; for (i=2;i(QS+1);i+) nextElement=arrayi; array0=nextElement; j=i; while (nextElementarrayj-1) arrayj=arrayj-1; j-; arrayj=nextElement; for (i=0;iQS;i+) Di=arrayQS-i; /

21、离的最近的四个数存储在Tem4数组中 double TemFix=0.0; for (i=0;iFix;i+) Temi=Di;/ cout 1/Temi endl; /提取速度V_EFix;V_NFix double V_EFix=0,V_NFix=0; for (i=0;iFix;i+) for (j=0;jQS;j+) if(Temi=D_BFj) V_Ei=Espeedj; V_Ni=Nspeedj; /计算E,N方向的速度 double Ev=0.0,Nv=0.0; double sum=0; for (i=0;iFix;i+) sum+=Temi; for (i=0;iFix;i+

22、) Ev+=Temi*V_Ei/sum; Nv+=Temi*V_Ni/sum; ouf1 setprecision(5) (100+n*0.5) (29+m*0.5) Ev Nv endl; Ev=0.0; Nv=0.0;/循环结束！ cout 计算完成 6。当数据点i ( xi , yi )到待定点P ( xP , yp )的距离di = sqr(xi2 - .yi2 ) i R时,该点即被选用。若选择的点数不够时,则应增大R的数值,直至数据点的个数n 满足要求。（3）选择二次曲面Z =Ax2+B xy +Cy2 +Dx + Ey +F 作为拟合面,则对应点的误差方程为vi = Ax2 +B

23、 xy +Cy2 +Dx + Ey +F - Zi 由n个数据点列出的误差方程为v = MX - Z X = (MTPM) - 1MTPZ 由于坐标原点移至该DEM格网点P ( xP , yp ) , 所以系数F就是待定点的内插值ZP 2.3 最小二乘配置法拟合水平速度场 最小二乘配置：最小二乘配置模型将观测量分为固定效应和随机效应两部分，能较好地解决含有倾向性、随机性因素的问题，适宜于系统误差的影响。最小二乘配置的模型为：式中其中，L为 n1维观测向量；AX为确定性参数部分；BS为随机效应部分；A为ntn维确定性参数系数矩阵； B为sts维随机效应部分系数矩阵；X为tn1维确定性参数向量；S为ts1信号向量；为 n1维观测误差向量。在本文中，确定性部分采用二次多项式，则tn为6；随机效应部分的系数 B 为单位矩阵，则ts=n,为点位纬度，为点位经度。 在最小二乘配置解式