摄影测量实验报告(前方交汇后方交汇)

摄影测量学实验报告学院： 地信院 班级： 测绘0904班 老师： 邹峥嵘 姓名： 张文佳 学号： 2011年11 月11 日空间后方交会空间前方交会程序编程实验一 实验目的1、 要求掌握运用摄影测量中空间后方交会-空间前方交会求解地面点的空间位置的方法和原理。2、 学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的计算，完成所给像对中两张像片各自的六个外方位元素的求解和精度评定。3、 根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，利用计算机编程语言前方交会编程，求解其对应的地面点在摄影测量坐标系中的坐标，从而达到通过摄影测量量测地面地理数据的目的。二 实验仪器1、 计算机2、 MATLAB计算机编程软件三、实验数据实验数据实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。此四对坐标运用最小二乘法求解左右像片的外方位元素， 即完成了空间后方的过程。另外还 给出了 5 对地面点在左右像片中的像平面坐标和左右像片的内方位元素。实验数据如下：f=150.000mm，x0=0，y0=0点号左片右片地面摄影测量坐标xyxyXYZGCP116.01279.963-73.9378.7065083.2055852.099527.925GCP288.5681.134-5.25278.1845780.025906.365571.549GCP313.362-79.37-79.122-78.8795210.8794258.446461.81GCP482.24-80.027-9.887-80.0895909.2644314.283455.484151.75880.555-39.95378.463214.618-0.231-76.0060.036349.88-0.782-42.201-1.022486.14-1.346-7.706-2.112548.035-79.962-44.438-79.736四、程序设计流程图1、后方交会输入GCP 的像点坐标 x，y计算旋转矩阵 R计算像点在像空间坐标系中的近似值xj, yj, 并组成误差方程的常数项,L=x-xj计算误差方程的系数项组成系数矩阵 A 组成法方程式计算系数 AA 常数项 AL 并求解外方位元素计算c、w、k、Xs、Ys、Zs 改正后的值确定初始值 c=w=k=0，Xs0，Ys0，Zs0小于跳出循环，完成迭代计算，精度评定判断改正数是否小于限差?大于此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（长度改正数小于0.01m,角度改正数小于0.0003，相当于1的角度值）为止。在这个过程中采用迭代计算的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。2、前方交会输入所需计算点的像平面坐标 x, y根据后方交会所得的旋转矩阵 Ra, Rb 计算像点在左、右像空间辅助坐标系中的坐标 Xa Ya Za，Xb Yb Zb计算摄影基线的三个坐标分量 Bx By Bz计算个点在左右像片中的的投影系数 Na ,Nb计算地面所求点在地面摄影测量坐标系中的坐标 X Y Z计算完毕七、实验原理公式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.,0.08856,0.,0.08224;y=0.,0.,-0.07937,-0.; %设置初始值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 %统计迭代计算次数%完成迭代计算，检验改正是是否符合要求while 1%计算旋转矩阵系数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 for n=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 endend%精度评定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.,0.,0.04988,0.08614,0.;ya=0.,-0.,-0.,-0.,-0.;xb=-0.,-0.,-0.,-0.,-0.;yb=0.,0.,-0.,-0.,-0.; %由外方位线元素计算基线分量 Bx, By, BzBx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1; %由外方位角元素计算像空间辅助坐标 X1, Y1, Z1 , X2, Y2, Z2R1=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.015Fb=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 , N2for n=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);EndX=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.0951 0.9950 -0.02910.0030 0.0287 0.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+003Ys0 = 5.0001e+003Zs0 = 2.0000e+003（5）左片迭代次数W = 5（6）左片外方元素协因素阵Qxx = 1.0e+009 * 1.3812 0.0062 0.4778 -0.0009 -0.0000 0.0001 0.0062 2.3416 0.1160 -0.0000 -0.0011 -0.0003 0.4778 0.1160 0.2659 -0.0003 -0.0001 0.0000 -0.0009 -0.0000 -0.0003 0.0000 0.0000 -0.0000 -0.0000 -0.0011 -0.0001 0.0000 0.0000 0.0000 0.0001 -0.0003 0.0000 -0.0000 0.0000 0.0000（7）左片单位权中误差mo =3.1794e-006（8）左片外方元素改正数中误差2.右片计算结果（1）旋转矩阵R = 0.9938 -0.1108 -0.0134 0.1101 0.9929 -0.0461 0.0184 0.0325 0.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+003Ys0 = 5.0879e+003Zs0 = 2.0299e+003（5）右片迭代次数W = 6（6）右片外方元素协因素阵Qxx = 1.0e+009 * 1.5207 0.1224 -0.4192 -0.0009 -0.0001 0.0001 0.1224 2.6049 -0.0578 -0.0001 -0.0013 0.0003 -0.4192 -0.0578 0.2084 0.0003 0.0000 -0.0000 -0.0009 -0.0001 0.0003 0.0000 0.0000 -0.0000 -0.0001 -0.0013 0.0000 0.0000 0.0000 -0.0000 0.0001 0.0003 -0.0000 -0.0000 -0.0000 0.0000（7）右片单位权中误差mo = 6.7173e-006（8）右片改正数中误差3.前方交会反算的地面摄影测量坐标X = 1.0e+003 * 5.3846 5.1324 5.4571 5.8154 5.5274Y = 1.0e+003 * 5.7446 5.0165 5.0414 5.0679 4.2922Z = 1.0e+003 *1.8901 1.8638 1.8633 1.8597 1.8368十、实验心得体会经过几个星期的努力，这个当初看似不可能完成的任务终于完成了，我感到很欣慰很有成就感。通过这次编程实验，加强了我对空间后方交会求解外方元素以及应用空间前方交会求解物点在摄影测量做坐标系中的三维坐标，同时让我认识到将理论联系实际的重要性。在本次试验中，我都遇到了很多困难，但都被我逐级突破。在这之前我并没