摄影测量空间前交后交资料.docx

摄影测量空间前交后交资料.docx

《摄影测量空间前交后交资料.docx》由会员分享，可在线阅读，更多相关《摄影测量空间前交后交资料.docx（16页珍藏版）》请在冰豆网上搜索。

摄影测量空间前交后交资料

空间后交-前交程序设计

（实验报告）

姓名：

班级：

学号：

时间：

空间后交-前交程序设计

一、实验目的

用C、VB或MATLAB语言编写空间后方交会-空间前方交会程序

⑴提交实习报告：

程序框图、程序源代码、计算结果、体会

⑵计算结果：

像点坐标、地面坐标、单位权中误差、外方位元素及其精度

二、实验数据

f=150.000mm，x0=0，y0=0

三、实验思路

1.利用空间后方交会求左右像片的外方位元素

（1）.获取m（于像片中选取两点，于地面摄影测量坐标系中选取同点，分别计算距离，距离比值即为m），x,y,f,X,Y,Z

（2）.确定未知数初始值Xs,Ys,Zs,q,w,k

（3）.计算旋转矩阵R

（4）.逐点计算像点坐标的近似值（x），（y）

（5）.组成误差方程式

（6）.组成法方程式

（7）.解求外方位元素

（8）.检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤（3）-（7）

2.利用求出的外方位元素进行空间前交，求出待定点地面坐标

（1）.用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1,R2

（2）.根据左、右像片的外方位元素，计算摄影基线分量Bx，By，Bz

（3）.计算像点的像空间辅助坐标（X1,Y1,Z1）和（X2,Y2,Z2）

（4）.计算点投影系数N1和N2

（5）.计算未知点的地面摄影测量坐标

四、实验过程

⑴程序框图

⑵程序代码

函数AandL

%求间接平差时需要的系数

%%%已知

%a=像点坐标x,b=像点坐标y，f内方位元素主距

%φ=q,ψ=w,κ=k

%像空间坐标系X,Y,Z

%地面摄影测量坐标系Xs,Ys,Zs

function[A1,L1,A2,L2]=AandL（a,b,f,q,w,k,X,Y,Z,Xs,Ys,Zs）

%%%%%%%%%%%选择矩阵元素

a1=cos（q）*cos（k）-sin（q）*sin（w）*sin（k）;

a2=-cos（q）*sin（k）-sin（q）*sin（w）*cos（k）;

a3=-sin（q）*cos（w）;

b1=cos（w）*sin（k）;

b2=cos（w）*cos（k）;

b3=-sin（w）;

c1=sin（q）*cos（k）+cos（q）*sin（w）*sin（k）;

c2=-sin（q）*sin（k）+cos（q）*sin（w）*cos（k）;

c3=cos（q）*cos（w）;

%%%%%%%共线方程的分子分母

X_=a1*（X-Xs）+b1*（Y-Ys）+c1*（Z-Zs）;

Y_=a2*（X-Xs）+b2*（Y-Ys）+c2*（Z-Zs）;

Z_=a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs）;

%%%%%%%近似值

x=-f*X_/Z_;

y=-f*Y_/Z_;

%%%%%%%A组成L组成

a11=1/Z_*（a1*f+a3*x）;

a12=1/Z_*（b1*f+b3*x）;

a13=1/Z_*（c1*f+c3*x）;

a21=1/Z_*（a2*f+a3*y）;

a22=1/Z_*（b2*f+b3*y）;

a23=1/Z_*（c2*f+c3*y）;

a14=y*sin（w）-（x/f*（x*cos（k）-y*sin（k））+f*cos（k））*cos（w）;

a15=-f*sin（k）-x/f*（x*sin（k）+y*cos（k））;

a16=y;

a24=-x*sin（w）-（y/f*（x*cos（k）-y*sin（k））-f*sin（k））*cos（w）;

a25=-f*cos（k）-y/f*（x*sin（k）+y*cos（k））;

a26=-x;

lx=a-x;

ly=b-y;

%%%%%%%%%组成一个矩阵，并返回

A1=[a11,a12,a13,a14,a15,a16];

A2=[a21,a22,a23,a24,a25,a26];

L1=lx;

L2=ly;

函数deg2dms

%%%%%%%%角度转度分秒

functiony=deg2dms（x）

a=floor（x）;

b=floor（（x-a）*60）;

c=（x-a-b/60）*3600;

y=a+（b/100）+（c/10000）;

函数dms2deg

%%%%%度分秒转度

functiony=dms2deg（x）

a=floor（x）;

b=floor（（x-a）*100）;

c=（x-a-b/100）*10000;

y=a+b/60+c/3600;

函数ok

%%%%%%%%%%%

%%%目的是为了保证各取的值的有效值

%%xy为n*1,a为1*n

functionresult=ok（xy,a）

formatshortg

i=size（xy,1）;

forn=1:

i

o=xy（n）-floor（xy（n,1））;

o=round（o*（10^a（n）））/（10^a（n））;

xy（n,1）=floor（xy（n,1））+o;

end

formatlongg

result=xy;

函数rad2dmsxy

%%%%求度分秒表现形式的三个外方位元素，三个角度

functionxydms=rad2dmsxy（xy）

[a,b,c,d,e,f]=testvar（xy）;

d=deg2dms（rad2deg（d））;

e=deg2dms（rad2deg（e））;

f=deg2dms（rad2deg（f））;

xydms=[a,b,c,d,e,f]';

函数spacehoujiao

%%%%%%%空间后交

%%%f

%%输入p（2*n,1）

%%像点坐标x,y,X,Y,Z,均为（n,1）

function[xy,m,R]=spacehoujiao（p,x,y,f,X,Y,Z）

formatlong;

%%%%%权的矢量化,这是等精度时的，如果非，将函数参数改为P

P=diag（p）;

%%求n

j=size（X,2）;

%%初始化

Xs=0;Ys=0;Zs=0;

forn=1:

j

Xs=Xs+X（n）;

Ys=Ys+Y（n）;

Zs=Zs+Z（n）;

end

Sx=sqrt（（x

（2）-x

（1））^2+（y

（2）-y

（1））^2）;%%%%两像点之间距离

Sd=sqrt（（X

（2）-X

（1））^2+（Y

（2）-Y

（1））^2）;%%%%两地面控制点之间距离

m=Sd/Sx;%%%%图像比例系数

Xs=Xs/j;

Ys=Ys/j;

Zs=m*f+Zs/j;

m0=0;

q=0;w=0;k=0;i=0;

a=rand（2*j,6）;

l=rand（2*j,1）;

%%%%

forn=1:

j

[a（2*n-1,:

）,l（2*n-1,1）,a（2*n,:

）,l（2*n,1）]=AandL（x（n）,y（n）,f,q,w,k,X（n）,Y（n）,Z（n）,Xs,Ys,Zs）;

end

det=inv（a'*P*a）*transpose（a）*P*l;

%%%%%%%%%循环体

while1

%%%%%%%%%%%%%%%%

[dXs,dYs,dZs,dq,dw,dk]=testvar（det）;

detXs=abs（dXs）;

detYs=abs（dYs）;

detZs=abs（dZs）;

detq=abs（dq）;

detw=abs（dw）;

detk=abs（dk）;

%%%%%%%%%

if（（detXs<0.01）&&（detYs<0.01）&&（detZs<0.01）&&（detq

break;

else

V=（a*det-l）;Q=inv（a'*P*a）;

m0=m0+sqrt（（V'*P*V）/（2*j-6））;%%m0需要每次的改正数算出来相加

%%%

Xs=Xs+dXs;

Ys=Ys+dYs;

Zs=Zs+dZs;

q=q+dq;

w=w+dw;

k=k+dk;

%%%

forn=1:

j

[a（2*n-1,:

）,l（2*n-1,1）,a（2*n,:

）,l（2*n,1）]=AandL（x（n）,y（n）,f,q,w,k,X（n）,Y（n）,Z（n）,Xs,Ys,Zs）;

end

det=inv（a'*P*a）*transpose（a）*P*l;

i=i+1;

%%%%

end

%%%

end

[dXs,dYs,dZs,dq,dw,dk]=testvar（det）;

detXs=abs（dXs）;

detYs=abs（dYs）;

detZs=abs（dZs）;

detq=abs（dq）;

detw=abs（dw）;

detk=abs（dk）;

V=（a*det-l）;

Q=inv（a'*P*a）;

m0=m0+sqrt（（V'*P*V）/（2*n-6））;

%%%

Xs=Xs+dXs;

Ys=Ys+dYs;

Zs=Zs+dZs;

q=q+dq;

w=w+dw;

k=k+dk;

%%%%%%%%%%%%%可以输出迭代次数的i

%%%%%%%%%%%%Xs,Ys,Zs,q,w,k,i,dXs,dYs,dZs,dq,dw,dk,detXs,detYs,detZs

%%%%%%%%%%%精度

mo=m0*sqrt（Q）;

m=[mo（1,1）,mo（2,2）,mo（3,3）,mo（4,4）,mo（5,5）,mo（6,6）]';

[mXs,mYs,mZs,mq,mw,mk]=testvar（m）;

%%%%%%%%%输出

xy=[Xs,Ys,Zs,q,w,k]';%%输出（6,1）的外方位元素

m=[m0,mXs,mYs,mZs,mq,mw,mk]';%%单位误差，各元素中误差

R=xyR（xy）;%%旋转矩阵

函数spaceqianjiao

%%空间前交

%输入f

%输入x1,y1,x2,y2,R1,R2,xy1,xy2（n,1）

%输出X,Y,Z（n,1）

function[X,Y,Z]=spaceqianjiao（x1,y1,x2,y2,f,R1,R2,xy1,xy2）

i=size（x1,2）;

[Xs1,Ys1,Zs1,q1,w1,k1]=testvar（xy1）;

[Xs2,Ys2,Zs2,q2,w2,k2]=testvar（xy2）;

forn=1:

i

[X1（n）,Y1（n）,Z1（n）]=testvar（R1*[x1（n）,y1（n）,-f]'）;

[X2（n）,Y2（n）,Z2（n）]=testvar（R2*[x2（n）,y2（n）,-f]'）;

Bx=Xs2-Xs1;

By=Ys2-Ys1;

Bz=Zs2-Zs1;

N1=（Bx*Z2（n）-Bz*X2（n））/（X1（n）*Z2（n）-X2（n）*Z1（n））;

N2=（Bx*Z1（n）-Bz*X1（n））/（X1（n）*Z2（n）-X2（n）*Z1（n））;

X（n）=Xs1+N1*X1（n）;

Z（n）=Zs1+N1*Z1（n）;

Y（n）=0.5*（（Ys1+N1*Y1（n））+（Ys2+N2*Y2（n）））;

end

函数testvar

%分割矩阵。

。

%%将矩阵的每行元素打包给元素。

。

%%%用法[Xs1,Ys1,Zs1,q1,w1,k1]=testvar（xy1）;

function[varargout]=testvar（arrayin）

fork=1:

nargout

varargout{k}=arrayin（k,:

）;

end

函数xyR

%计算旋转矩阵，通过六个内方位元素

%xy（6,1）

function[R]=xyR（xy）

[a,b,c,q,w,k]=testvar（xy）;

a1=cos（q）*cos（k）-sin（q）*sin（w）*sin（k）;

a2=-cos（q）*sin（k）-sin（q）*sin（w）*cos（k）;

a3=-sin（q）*cos（w）;

b1=cos（w）*sin（k）;

b2=cos（w）*cos（k）;

b3=-sin（w）;

c1=sin（q）*cos（k）+cos（q）*sin（w）*sin（k）;

c2=-sin（q）*sin（k）+cos（q）*sin（w）*cos（k）;

c3=cos（q）*cos（w）;

R=[a1,a2,a3;b1,b2,b3;c1,c2,c3];

命令行

f=0.15;

%%空间后交

%像片坐标

x11=（1e-3）*[16.012,88.56,13.362,82.24];

y11=（1e-3）*[79.963,81.134,-79.37,-80.027];

x21=（1e-3）*[-73.93,-5.252,-79.122,-9.887];

y21=（1e-3）*[78.706,78.184,-78.879,-80.089];

%地面摄影测量坐标

X1=[5083.205,5780.02,5210.879,5909.264];

Y1=[5852.099,5906.365,4258.446,4314.283];

Z1=[527.925,571.549,461.81,455.484];

p=[1,1,1,1,1,1,1,1];

%%%

%xy1,xy2六个内方位元素矩阵，xy1,左片xy2，右片

%m1,m2误差矩阵，m1,左片m2，右片

%R1,R2旋转矩阵R1,左片R2，右片

[xy1,m1,R1]=spacehoujiao（p,x11,y11,f,X1,Y1,Z1）;

[xy2,m2,R2]=spacehoujiao（p,x21,y21,f,X1,Y1,Z1）;

xy1dms=ok（rad2dmsxy（xy1）,[3,3,3,4,4,4]）

m1

xy2dms=ok（rad2dmsxy（xy2）,[3,3,3,4,4,4]）

m2

%%%%空间前交

x12=（1e-3）*[51.758,14.618,49.88,86.14,48.035];

y12=（1e-3）*[80.555,-0.231,-0.782,-1.346,-79.962];

x22=（1e-3）*[-39.953,-76.006,-42.201,-7.706,-44.438];

y22=（1e-3）*[78.463,0.036,-1.022,-2.112,-79.736];

[X2,Y2,Z2]=spaceqianjiao（x12,y12,x22,y22,f,R1,R2,xy1,xy2）;

X2d=（ok（X2',[3,3,3,3,3]））'

Y2d=（ok（Y2',[3,3,3,3,3]））'

Z2d=（ok（Z2',[3,3,3,3,3]））'

（3）计算结果

左片六个外方位元素xy1dms单位中误差和六个外方位元素中误差m1

左片六个外方位元素xy2dms单位中误差和六个外方位元素中误差m2

空间前交得到的坐标X2d，Y2d,Z2d

五、心得体会

这次空间后交-前交程序设计，我用matlab语言进行编程。

在矩阵的计算方面，matlab较c++优势明显，但是，由于自己matlab懂得不多，同时，不明白matlab的调试，导致了我不得不花费很长的时间寻找打错的符号、公式。

此外，在进行空间前方交会的编程时，对循环的判定，循环体中的公式，中误差的计算，都不清楚，这无疑降低了编程的速度。

在这次程序设计中，出现了非常多的问题，比如，对三个角外方位元素进行转度分秒时，由于前、后方交会的程序相对独立，没有考虑到对后方造成的影响，会使结果偏差特别大。

花费了大量时间后，我终于将前、后方交会的程序编出来了。

这很令人喜悦，但在某种程度上，也是对自己的一种鞭策。

花费了这么多时间编出来的程序，真是渣的不行；在编程时，效率十分低下，看看视频，听听音乐，一个小时就过去了，代码才敲了几行，还有错误；写了这么久的空间前、后交会，只知道大致的原理，具体一点的就一问三不知了。

总之，在这次试验中，自己虽然牺牲了很多的课余时间，但是收获还是非常多的。

同时，自己也发现，编程能力的提高，在明白基本的知识后，需要通过大量的练习进行提高。