摄影测量空间前交后交资料.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）&&（detqbreak;
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的调试，导致了我不得不花费很长的时间寻找打错的符号、公式。
此外，在进行空间前方交会的编程时，对循环的判定，循环体中的公式，中误差的计算，都不清楚，这无疑降低了编程的速度。
在这次程序设计中，出现了非常多的问题，比如，对三个角外方位元素进行转度分秒时，由于前、后方交会的程序相对独立，没有考虑到对后方造成的影响，会使结果偏差特别大。
花费了大量时间后，我终于将前、后方交会的程序编出来了。
这很令人喜悦，但在某种程度上，也是对自己的一种鞭策。
花费了这么多时间编出来的程序，真是渣的不行；在编程时，效率十分低下，看看视频，听听音乐，一个小时就过去了，代码才敲了几行，还有错误；写了这么久的空间前、后交会，只知道大致的原理，具体一点的就一问三不知了。
总之，在这次试验中，自己虽然牺牲了很多的课余时间，但是收获还是非常多的。
同时，自己也发现，编程能力的提高，在明白基本的知识后，需要通过大量的练习进行提高。
升级会员 