摄影测量实验报告前方交汇后方交汇.doc
- 文档编号:1307204
- 上传时间:2023-04-30
- 格式:DOC
- 页数:12
- 大小:168KB
摄影测量实验报告前方交汇后方交汇.doc
《摄影测量实验报告前方交汇后方交汇.doc》由会员分享,可在线阅读,更多相关《摄影测量实验报告前方交汇后方交汇.doc(12页珍藏版)》请在冰点文库上搜索。
摄影测量学
实验报告
学院:
地信院
班级:
测绘0904班
老师:
邹峥嵘
姓名:
张文佳
学号:
0405090921
2011年11月11日
空间后方交会——空间前方交会
程序编程实验
一.实验目的
1、要求掌握运用摄影测量中空间后方交会-空间前方交会求解地面点的空间位置的方法和原理。
2、学会运用空间后方交会的原理,根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标,利用计算机编程语言实现空间后方交会的计算,完成所给像对中两张像片各自的六个外方位元素的求解和精度评定。
3、根据空间后方交会所得的两张像片的内外方位元素,利用同名像点在左右像片上的坐标,利用计算机编程语言前方交会编程,求解其对应的地面点在摄影测量坐标系中的坐标,从而达到通过摄影测量量测地面地理数据的目的。
二.实验仪器
1、计算机
2、MATLAB计算机编程软件
三、实验数据
实验数据实验数据包含四个地面控制点(GCP)的地面摄影测量坐标及在左右像片中的像平面坐标。
此四对坐标运用最小二乘法求解左右像片的外方位元素,即完成了空间后方的过程。
另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。
实验数据如下:
f=150.000mm,x0=0,y0=0
点号
左片
右片
地面摄影测量坐标
x
y
x
y
X
Y
Z
GCP1
16.012
79.963
-73.93
78.706
5083.205
5852.099
527.925
GCP2
88.56
81.134
-5.252
78.184
5780.02
5906.365
571.549
GCP3
13.362
-79.37
-79.122
-78.879
5210.879
4258.446
461.81
GCP4
82.24
-80.027
-9.887
-80.089
5909.264
4314.283
455.484
1
51.758
80.555
-39.953
78.463
2
14.618
-0.231
-76.006
0.036
3
49.88
-0.782
-42.201
-1.022
4
86.14
-1.346
-7.706
-2.112
5
48.035
-79.962
-44.438
-79.736
四、程序设计流程图
1、后方交会
输入GCP的像点坐标x,y
计算旋转矩阵R
计算像点在像空间坐标系中的近似值xj,yj,并组成误差方程的常数项,L=x-xj
计算误差方程的系数项组成系数矩阵A
组成法方程式
计算系数A’A常数项A’L
并求解外方位元素
计算c、w、k、Xs、Ys、Zs改正后的值
确定初始值c=w=k=0,Xs0,Ys0,Zs0
小于
跳出循环,完成迭代计算,精度评定
判断改正数是否小于限差?
大于
此过程完成空间后方交会求解像片的外方位元素,其中改正数小于限差(长度改正数小于0.01m,角度改正数小于0.0003,相当于1’的角度值)为止。
在这个过程中采用迭代计算的方法,是外方位元素逐渐收敛于理论值,每次迭代所得的改正数都应加到上一次的初始值之中。
2、前方交会
输入所需计算点的像平面坐标x,y
根据后方交会所得的旋转矩阵Ra,Rb计算像点在左、右像空间辅助坐标系中的坐标XaYaZa,XbYbZb
计算摄影基线的三个坐标分量BxByBz
计算个点在左右像片中的的投影系数Na,Nb
计算地面所求点在地面摄影测量坐标系中的坐标XYZ
计算完毕
七、实验原理公式
1、后方交会中运用的共线方程数学模型
3、前方交会与后方交会中均用到旋转矩阵进行的坐标转换
4、精度评定中均采用最小二乘准则进行平差计算
5.前方交会的转换
八.实验源程序
1.空间后方交会(以左片为例)
%已知地面摄影测量坐标
Xg=[5083.205,5780.02,5210.879,5909.264];
Yg=[5852.099,5906.365,4258.466,4314.283];
Zg=[527.925,571.549,461.81,455.484];
%对应地面坐标的左片像点坐标
x=[0.016012,0.08856,0.013362,0.08224];
y=[0.079963,0.081134,-0.07937,-0.080027];
%设置初始值
c=0;w=0;k=0;
f=0.15;
Dg=sqrt((Xg
(1)-Xg
(2))^2+(Yg
(1)-Yg
(2))^2);
Ds=sqrt((x
(1)-x
(2))^2+(y
(1)-y
(2))^2);
p=Dg/Ds;
Xs0=1/4*(Xg
(1)+Xg
(2)+Xg(3)+Xg(4));
Ys0=1/4*(Yg
(1)+Yg
(2)+Yg(3)+Yg(4));
Zs0=p*f+(Zg
(1)+Zg
(2)+Zg(3)+Zg(4))/4;
W=0%统计迭代计算次数
%完成迭代计算,检验改正是是否符合要求
while1
%计算旋转矩阵系数
a1=cos(c)*cos(k)-sin(c)*sin(w)*sin(k);
a2=-cos(c)*sin(k)-sin(c)*sin(w)*cos(k);
a3=-sin(c)*cos(w);
b1=cos(w)*sin(k);
b2=cos(w)*cos(k);
b3=-sin(w);
c1=sin(c)*cos(k)+cos(c)*sin(w)*sin(k);
c2=-sin(c)*cos(k)+cos(c)*sin(w)*cos(k);
c3=cos(c)*cos(w);
R=[a1,a2,a3;b1,b2,b3;c1,c2,c3]
forn=1:
1:
4
X(n)=a1*(Xg(n)-Xs0)+b1*(Yg(n)-Ys0)+c1*(Zg(n)-Zs0);
Y(n)=a2*(Xg(n)-Xs0)+b2*(Yg(n)-Ys0)+c2*(Zg(n)-Zs0);
Z(n)=a3*(Xg(n)-Xs0)+b3*(Yg(n)-Ys0)+c3*(Zg(n)-Zs0);
xj(n)=-f*X(n)/Z(n);%计算近似点坐标
yj(n)=-f*Y(n)/Z(n);
a11(n)=1/Z(n)*(a1*f+a3*x(n));%矩阵系数
a12(n)=1/Z(n)*(b1*f+b3*x(n));
a13(n)=1/Z(n)*(c1*f+c3*x(n));
a21(n)=1/Z(n)*(a2*f+a3*y(n));
a22(n)=1/Z(n)*(b2*f+b3*y(n));
a23(n)=1/Z(n)*(c2*f+c3*y(n));
a14(n)=y(n)*sin(w)-(x(n)/f*(x(n)*cos(k)-y(n)*sin(k))+f*cos(k))*cos(w);
a15(n)=-f*sin(k)-x(n)/f*(x(n)*sin(k)+y(n)*cos(k));
a16(n)=y(n);
a24(n)=-x(n)*sin(w)-(y(n)/f*(x(n)*cos(k)-y(n)*sin(k))-f*sin(k))*cos(w);
a25(n)=-f*cos(k)-y(n)/f*(x(n)*sin(k)+y(n)*cos(k));
a26(n)=-x(n);
end
A=[a11
(1),a12
(1),a13
(1),a14
(1),a15
(1),a16
(1);
a21
(1),a22
(1),a23
(1),a24
(1),a25
(1),a26
(1);
a11
(2),a12
(2),a13
(2),a14
(2),a15
(2),a16
(2);
a21
(2),a22
(2),a23
(2),a24
(2),a25
(2),a26
(2);
a11(3),a12(3),a13(3),a14(3),a15(3),a16(3);
a21(3),a22(3),a23(3),a24(3),a25(3),a26(3);
a11(4),a12(4),a13(4),a14(4),a15(4),a16(4);
a21(4),a22(4),a23(4),a24(4),a25(4),a26(4)];
L=[x
(1)-xj
(1);y
(1)-yj
(1);x
(2)-xj
(2);y
(2)-yj
(2);x(3)-xj(3);y(3)-yj(3);x(4)-xj(4);y(4)-yj(4)];
Xp=inv((A')*A)*(A')*L
V=A*Xp-L
Xs0=Xs0+Xp(1,1)
Ys0=Ys0+Xp(2,1)
Zs0=Zs0+Xp(3,1)
c=c+Xp(4,1);w=w+Xp(5,1);k=k+Xp(6,1);
W=W+1
%判断收敛条件if((abs(Xp(1,1))<0.01&&abs(Xp(2,1))<0.01&&abs(Xp(3,1))<0.01&&abs(Xp(4,1))<0.0003&&abs(Xp(5,1))<0.0003&&abs(Xp(6,1))<0.0003))
break
end
end
%精度评定
Qxx=inv((A')*A)%外方元素协因素阵
mo=sqrt(((V)'*V)/(2*8-6))%单位权中误差
mi=mo*sqrt(Qxx)%外方元素改正数中误差
2.前方交会
%后方交会中计算出的改正后外方元素
Xs1=4999.7;Ys1=5000.1;Zs1=2000.0;
Xs2=5896.3;Ys2=5087.9;Zs2=2029.9;
%量测像点左右片坐标
xa=[0.051758,0.014618,0.04988,0.08614,0.048035];
ya=[0.080555,-0.000231,-0.000782,-0.001346,-0.079962];
xb=[-0.039953,-0.076006,-0.042201,-0.007706,-0.044438];
yb=[0.078463,0.000036,-0.001022,-0.002112,-0.079736];
%由外方位线元素计算基线分量Bx,By,Bz
Bx=Xs2-Xs1;
By=Ys2-Ys1;
Bz=Zs2-Zs1;
%由外方位角元素计算像空间辅助坐标X1,Y1,Z1,X2,Y2,Z2
R1=[0.9955,-0.0951,-0.0002;0.0951,0.9950,-0.0291;0.0030,0.0287,0.9996];
R2=[0.9938,-0.1108,-0.0134;0.1101,0.9929,-0.0461;0.0184,0.0325,0.9988];
Fa=R1*[xa
(1),xa
(2),xa(3),xa(4),xa(5);ya
(1),ya
(2),ya(3),ya(4),ya(5);-0.015,-0.015,-0.015,-0.015,-0.015]
Fb=R2*[xb
(1),xb
(2),xb(3),xb(4),xb(5);yb
(1),yb
(2),yb(3),yb(4),yb(5);-0.015,-0.015,-0.015,-0.015,-0.015]
Xa=[Fa(1,1),Fa(1,2),Fa(1,3),Fa(1,4),Fa(1,5)];
Ya=[Fa(2,1),Fa(2,2),Fa(2,3),Fa(2,4),Fa(2,5)];
Za=[Fa(3,1),Fa(3,2),Fa(3,3),Fa(3,4),Fa(3,5)];
Xb=[Fb(1,1),Fb(1,2),Fb(1,3),Fb(1,4),Fb(1,5)];
Yb=[Fb(2,1),Fb(2,2),Fb(2,3),Fb(2,4),Fb(2,5)];
Zb=[Fb(3,1),Fb(3,2),Fb(3,3),Fb(3,4),Fb(3,5)];
%计算点投影系数N1,N2
forn=1:
1:
5
Na(n)=(Bx*Zb(1,(n))-Bz*Xb(1,(n)))/(Xa(1,(n))*Zb(1,(n))-Xb(1,(n))*Za(1,(n)));
Nb(n)=(Bx*Zb(1,(n))-Bz*Xb(1,(n)))/(Xa(1,(n))*Zb(1,(n))-Xb(1,(n))*Za(1,(n)));
%计算地面坐标XA,YA,ZA
Xg(n)=Na(n)*Xa(n)+Xs1;
Zg(n)=Na(n)*Za(n)+Zs1;
Yg(n)=1/2*((Ys1+Na(n)*Ya(n))+(Ys2+Nb(n)*Yb(n)));
End
X=[Xg
(1),Xg
(2),Xg(3),Xg(4),Xg(5)]
Y=[Yg
(1),Yg
(2),Yg(3),Yg(4),Yg(5)]
Z=[Zg
(1),Zg
(2),Zg(3),Zg(4),Zg(5)]
九、计算结果,精度评定
1.左片计算结果
(1)旋转矩阵
R=0.9955-0.0951-0.0002
0.09510.9950-0.0291
0.00300.02870.9996
(2)左片外方元素改正数
Xp=
1.0e-003*
0.0277,-0.4099,0.0058,0.0000,0.0000,0.0000
(3)左片像点坐标改正数
V=
1.0e-005*
-0.2514,0.3595,-0.0704,-0.2227,0.5448,0.3856,-0.2179,-0.5205
(4)左片外方元素改正结果
Xs0=4.9997e+003
Ys0=5.0001e+003
Zs0=2.0000e+003
(5)左片迭代次数
W=5
(6)左片外方元素协因素阵
Qxx=
1.0e+009*
1.38120.00620.4778-0.0009-0.00000.0001
0.00622.34160.1160-0.0000-0.0011-0.0003
0.47780.11600.2659-0.0003-0.00010.0000
-0.0009-0.0000-0.00030.00000.0000-0.0000
-0.0000-0.0011-0.00010.00000.00000.0000
0.0001-0.00030.0000-0.00000.00000.0000
(7)左片单位权中误差
mo=3.1794e-006
(8)左片外方元素改正数中误差
2.右片计算结果
(1)旋转矩阵
R=0.9938-0.1108-0.0134
0.11010.9929-0.0461
0.01840.03250.9988
(2)右片外方元素改正数
Xp=
-0.0002,0.0072,-0.0003,0.0000,0.0000,0.0000
(3)右片像点坐标改正数
V=
1.0e-004*
-0.0639,-0.0103,0.0157,-0.0482,-0.0735,0.1193,0.1226,-0.0607
(4)右片外方元素改正结果
Xs0=5.8963e+003
Ys0=5.0879e+003
Zs0=2.0299e+003
(5)右片迭代次数
W=6
(6)右片外方元素协因素阵
Qxx=
1.0e+009*
1.52070.1224-0.4192-0.0009-0.00010.0001
0.12242.6049-0.0578-0.0001-0.00130.0003
-0.4192-0.05780.20840.00030.0000-0.0000
-0.0009-0.00010.00030.00000.0000-0.0000
-0.0001-0.00130.00000.00000.0000-0.0000
0.00010.0003-0.0000-0.0000-0.00000.0000
(7)右片单位权中误差
mo=6.7173e-006
(8)右片改正数中误差
3.前方交会
反算的地面摄影测量坐标
X=
1.0e+003*
5.38465.13245.45715.81545.5274
Y=
1.0e+003*
5.74465.01655.04145.06794.2922
Z=
1.0e+003*
1.89011.86381.86331.85971.8368
十、实验心得体会
经过几个星期的努力,这个当初看似不可能完成的任务终于完成了,我感到很欣慰很有成就感。
通过这次编程实验,加强了我对空间后方交会求解外方元素以及应用空间前方交会求解物点在摄影测量做坐标系中的三维坐标,同时让我认识到将理论联系实际的重要性。
在本次试验中,我都遇到了很多困难,但都被我逐级突破。
在这之前我并没有系统性的学习过MATLAB,甚至可以说是根本没用过,虽然通过了计算机二级,但如果仅凭自己的能力,用C编写这个程序,还是觉得毫无头绪。
后来听同学说MATLAB简单易学,容易掌握,而且处理矩阵运算有强大的优势,所以我毅然选择了MATLAB作为本次实验编程的软件。
花了几天的时间学习了MATLAB的编绘语言以及运算法则。
刚开始编写时毫无头绪,而且觉得变量多得错综复杂,搞得我两眼发晕。
后来我没有急着去编写,而是回归到课本,将空间后方交会和前方交会有关的理论知识系统认真的看了几遍,了解了整个计算过程后,弄清变量之间的转换关系后,再将整个计算过程绘制成流程图,再进行编程。
完成了有关代码的基本编写后,就开始了程序的调试了,这是让我最为痛苦的一个环节,看似简单,而我却花了很长时间,语法错误比较容易发现,根据软件提示是很容易找出来的,其中很容易出现的错误是错误的拼写,矩阵维数不一致等问题。
就在没有语法错误的时候我满怀希望的再次运行,却得到的了无限循环,让我很吃惊,强制退出程序后,重新检查了很多遍源程序,把书上的公式和自己的进行了详细比对,最后才发现因为马虎导致有的公式漏打括号,变量输错等原因,进行改正后得到了正确结果。
在整个过程中,因为时间等原因,我没有对MATLAB进行深入学习,没有运用到程序函数功能,用了两个m文件进行了左片右片的程序编写,导致了很多重复步骤,效率低下,也没有注意到程序运行结果的美观和清晰明了程度。
在下次进行运用时会多加注意。
总之,整个实验算是完成了,我不仅对后方交会和前方交会的实用原理更加明了,初步学习了MATLAB,通过本次实验也增强了我分析问题和处理问题的能力,也给了我一个机会把理论知识运用到实际计算中去的机会。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 摄影 测量 实验 报告 前方 交汇 后方