速度场模型.docx
《速度场模型.docx》由会员分享，可在线阅读，更多相关《速度场模型.docx（20页珍藏版）》请在冰豆网上搜索。
速度场模型
空间大地测量学作业
学院：
测绘科学与技术学院
专业：
测绘工程
姓名：
曹延虎
学号：
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个数据，如下图红线所示。
图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:
how2
forj2=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）];
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）;
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局部多项式法拟合水平速度场
多项式的形式有三种：
一次多项式、二次多项式和三次多项式。
在插值时，寻找一个适合的函数并非易事，如多项式阶数太大的问题，其波动也会很大。
针对这个问题，局部多项式法是一个不错的选择。
它作为常用的拟合方法之一，能在给定搜索区域内对插值对象所有点插值出合理特定阶数的多项式，局部多项式插值产生的曲面根多依赖于局部的变异。
常见的多项式形式如下：
本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。
图4拟合效果图
1.2.1Matlab源程序（原始数据及中间值见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）;
fori3=1:
how3
forj3=1:
10
ifj3==1
array3（i3,j3）=1;%1
elseifj3==2
array3（i3,j3）=X（i3）;%X
elseifj3==3
array3（i3,j3）=Y（i3）;%Y
elseifj3==4
array3（i3,j3）=X（i3）*Y（i3）;%X*Y
elseifj3==5
array3（i3,j3）=X（i3）^2;%X^2
elseifj3==6
array3（i3,j3）=Y（i3）^2;%Y^2
elseifj3==7
array3（i3,j3）=X（i3）^2*Y（i3）;%X*Y^2
elseifj3==8
array3（i3,j3）=Y（i3）^2*X（i3）;%Y*X^2
elseifj3==9
array3（i3,j3）=X（i3）^3;%X^3
elseifj3==10
array3（i3,j3）=Y（i3）^3;%Y^3
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
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）;
%%计算参数的中误差
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'）;%读入待求点的经度Y
Unl3_X=xlsread（'原始数据经纬度.xls',1,'E127:
E157'）;%读入待求点的纬度X
Oritation=zeros（size（Unb3_Y）,10）;
fork=1:
size（Unb3_Y）
forn=1:
10
ifn==1
Oritation（k,n）=1;%1
elseifn==2
Oritation（k,n）=Unl3_X（k）;%X
elseifn==3
Oritation（k,n）=Unb3_Y（k）;%Y
elseifn==4
Oritation（k,n）=Unl3_X（k）*Unb3_Y（k）;%X*Y
elseifn==5
Oritation（k,n）=Unl3_X（k）^2;%X^2
elseifn==6
Oritation（k,n）=Unb3_Y（k）^2;%Y^2
elseifn==7
Oritation（k,n）=Unl3_X（k）^2*Unb3_Y（k）;%X*Y^2
elseifn==8
Oritation（k,n）=Unb3_Y（k）^2*Unl3_X（k）;%Y*X^2
elseifn==9
Oritation（k,n）=Unl3_X（k）^3;%X^3
elseifn==10
Oritation（k,n）=Unb3_Y（k）^3;%Y^3
end
end
end
X_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点的插值。
其中
，表示由离散点到待插值点的距离。
计算格网点数据时，给每个观测点赋予一个权重，所有点的权重之和为1.0。
.如果某个观测点与格网点重合，那么这个观测点的权重为1.0，其它观测点的权重均为0.0。
该方法的一个特点就是在格网区域内，观测点周围将为产生“牛眼”现象。
我们可以通过一定的方法进行消除，比如在进行格网化的时候，可以选择一个圆滑参数，该参数只要大于零，则能保证对于一个特定的节点，所有的观测点都不能被赋予权重1.0，包括该点与结点重合。
本次拟合水平速度场所采用的数据与多面函数法的已知数据一样，本次实验也拟合了30个数据，如下图红线所示。
图5拟合效果图
1.3.1C++源程序（原始数据见文件夹）
#include
#include
#include
#include
usingnamespacestd;
#defineinfile1"C:
\\Users\\pc\\Desktop\\反距离加权法\\原始数据.txt"//起算数据文件
#defineoutfile1"C:
\\Users\\pc\\Desktop\\反距离加权法\\拟合数据.txt"//未知点速度场文件
constintQS=40;//读取起算数据的个数
constintFix=5;//需要参考点个数
intmain（）
{
ifstreaminf1;
ofstreamouf1;
inf1.open（infile1）;
if（inf1.fail（））
cout<<"cannotopeninfile1:
"<ouf1.open（outfile1）;
if（ouf1.fail（））
cout<<"cannotopenoutfile1:
"<//定义循环变量
inti=0,j=0,m=0,n=0;
//定义b,l,Espeed,Nspeed数组变量
floatl[QS],b[QS],Espeed[QS],Nspeed[QS];
//从文件读取起算点数据
inf1>>l[0]>>b[0]>>Espeed[0]>>Nspeed[0];
while（!
inf1.eof（））
{
++i;
inf1>>l[i]>>b[i]>>Espeed[i]>>Nspeed[i];
}
//定义未知点
doublelo=0,bo=0;
//计算待插值点，循环开始！
！
！
！
！
！
！
！
！
！
！
！
！
！
for（m=0;m<13;m++）
for（n=0;n<21;n++）
{
lo=100+0.5*n,bo=29+0.5*m;
//定义存储待插值点到起算点的距离数组
doubleD[QS]={0.0},D_BF[QS]={0.0};
//计算待插值点到起算点的距离
for（i=0;i{
D[i]=1/（（lo-l[i]）*（lo-l[i]）+（bo-b[i]）*（bo-b[i]））;
D_BF[i]=D[i];
}
//对数组D进行从大到小排序
doublenextElement=0;
doublearray[QS+1]={0.0};
for（i=0;i{
array[i+1]=D[i];
}
i=0,j=0;
for（i=2;i<（QS+1）;i++）
{
nextElement=array[i];
array[0]=nextElement;
j=i;
while（nextElement{
array[j]=array[j-1];
j--;
}
array[j]=nextElement;
}
for（i=0;i{
D[i]=array[QS-i];
}
//离的最近的四个数存储在Tem[4]数组中
doubleTem[Fix]={0.0};
for（i=0;i{
Tem[i]=D[i];
//cout<<1/Tem[i]<}
//提取速度V_E[Fix];V_N[Fix]
doubleV_E[Fix]={0},V_N[Fix]={0};
for（i=0;ifor（j=0;j{
if（Tem[i]==D_BF[j]）
{
V_E[i]=Espeed[j];
V_N[i]=Nspeed[j];
}
}
//计算E,N方向的速度
doubleEv=0.0,Nv=0.0;
doublesum=0;
for（i=0;i{
sum+=Tem[i];
}
for（i=0;i{
Ev+=Tem[i]*V_E[i]/sum;
Nv+=Tem[i]*V_N[i]/sum;
}
ouf1<Ev=0.0;
Nv=0.0;
}
//循环结束！
！
！
！
！
！
！
！
！
！
！
！
！
！
！
！
！
cout<<"计算完成"<inf1.close（）;
ouf1.close（）;
}
二其他方法简介
2.1欧拉矢量法拟合水平速度场
欧拉矢量法是最常用的一种速度场建模方法，但是需要假定建立模型的区域为刚体。
2.1块体欧拉矢量法
如果一个刚体绕固定点旋转,则该刚体由某一给定位置到另一位置的滑移,等效于绕通过该固定点的某一固定旋转轴旋转。
如果把地心当作固定点,块体的水平运动可以视为绕通过地心的定轴旋转。
假定划分的块体为刚性体,块体作一致性欧拉运动,内部站点间的几何构造在这个块体运动中保持不变。
假定划分的每个块体为刚性,提取各块体内的站点坐标及速度值,计算各块体的欧拉矢量,对站点速度进行插值。
计算方法如下：
设
为某块体上的站心坐标速度，
与该块体欧拉矢量式
为：
（2-1）
式中，r为地球平均半径,
为地理坐标经度，
为地理坐标纬度；
分别为站点东方向和北方向速度，
为欧拉矢量。
构造的误差方程为：
式中，A为欧拉矢量系数矩阵，X为欧拉矢量，L为建模点E和N两个方向的速度组成的向量，n为该局部区域内速度已知点的个数。
根据最小二乘原理得：
其中：
计算欧拉矢量后，利用式（2-1）对待插值点进行计算。
2.2局域欧拉矢量法
目前，对于中国大陆块体划分还没有统一的方法，在块体边界处一直存在争议，且无论怎么划分块体，都无法保证所划分块体为一个刚性体。
但有一点可以肯定，周部小的区域处在同一块体的可能性比较大，所以采用局部数据求欧拉矢量的方法可以提高欧拉矢量模型的可靠性。
具体做法如下：
1）在待插值点一定范围内搜索建模点，当搜索到的建模点数在4-10之间，利用搜索到的建模点求取欧拉矢量，并计算待插值点速度。
如果一定范围内建模点数大于10个，当搜索到10个建模点时，停止搜索，利用搜索到的10个建模点,求取欧拉矢量,并计算待插值点速度。
如果搜索到的建模点数少于4个，进行第二步。
2）在待插值点一定范围内搜索建模点，如果建模点数大于4个，选择距离待插值点最近的四个点求取欧拉矢量，并计算待插值点速度，如果建模点数少于4个,则此待插值点不适合采用局域欧拉矢量法计算其速度。
在利用欧拉矢量法建立待求点局部区域速度场模型时如何确定待求点周围局部区域点是使用该方法的关键。
2.2移动曲面拟合水平速度场
移动曲面拟合法是一种局部逼近方法，利用内插点周围的数据建立拟合曲面，使其到各数据点的距离之加权平方和为极小，而这个曲面在内插点上的值就是所求的内插值。
设Ｐ为要内插的点，对Ｐ构造相应的曲面：
为了确定曲面系数，需要选取Ｐ点周围的数据点。
选取的方法通常是以Ｐ为中心，Ｒ为半径的圆内的数据点，凡落在圆内的点即被选用，所选取的点数根据所采用的曲面函数来确定，以使得所选点的个数大于曲面方程系数的个数。
当点数不够时，则扩大Ｒ的值。
设所选取点的坐标为（ｘｉ，ｙｉ）（ｉ＝１，２，…，ｎ，ｎ≥６），且设Ｐ的坐标为（ｘｐ，ｙｐ）。
将（ｘｉ，ｙｉ）改化到以Ｐ为原点的局部坐标系中，即
则由Ｎ个数据点的值，可得
对每个数据点赋以权ωｉ，这里ωｉ不代表数据点的观测精度，而是反映该数据点与内插相关程度。
因此，对于权ωｉ确定的原则应与该数据点与内插点的距离有关，距离越小，它内插点的影响越大，则权应越大。
由最小二乘法解如下的带权极小问题：
由此解得系数ｃｉ，ｊ，从而得到Ｐ所对应的二次曲面方程，进而得到所求的内插点的值。
2、移动曲面拟合步骤：
（1）根据实际内插要求,解算待定点P的平面坐标（xP,yp）
（2）为了选取邻近的数据点,以待定点P为圆心,以R为半径作圆,凡落在圆内的数据点即被选用。
在二次曲面内插时,考虑到计算方便,将坐标原点移至该DEM格网点P（xP,yp）由于二次曲面系数个数为6,要求选用的数据点个数n>6。
当数据点i（xi,yi）到待定点P（xP,yp）的距离di=sqr（xi2-.yi2）i若选择的点数不够时,则应增大R的数值,直至数据点的个数n满足要求。
（3）选择二次曲面Z=Ax2+Bxy+Cy2+Dx+Ey+F作为拟合面,则对应点的误差方程为vi=Ax2+Bxy+Cy2+Dx+Ey+F-Zi由n个数据点列出的误差方程为v=MX-ZX=（MTPM）-1MTPZ由于坐标原点移至该DEM格网点P（xP,yp）,所以系数F就是待定点的内插值ZP
2.3最小二乘配置法拟合水平速度场
最小二乘配置：
最小二乘配置模型将观测量分为固定效应和随机效应两部分，能较好地解决含有倾向性、随机性因素的问题，适宜于系统误差的影响。
最小二乘配置的模型为：
式中
其中，L为n×1维观测向量；AX为确定性参数部分；BS为随机效应部分；A为n×tn维确定性参数系数矩阵；B为s×ts维随机效应部分系数矩阵；X为tn×1维确定性参数向量；S为ts×1信号向量；Δ为n×1维观测误差向量。
在本文中，确定性部分采用二次多项式，则tn为6；随机效应部分的系数B为单位矩阵，则ts=n,ψ为点位纬度，λ为点位经度。
在最小二乘配置解式
升级会员 