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

摄影测量学实验报告学院：地信院班级：测绘0904班老师：邹峥嵘姓名：张文佳学号：04050909212011年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、后方交会此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（长度改计算旋转矩阵R计算像点在像空间坐标系中的近似值xj,yj,并组成误差方程的常数项,L=x-xj计算误差方程的系数项组成系数矩阵A组成法方程式计算系数A’A常数项A’L并求解外方位元素计算c、w、k、Xs、Ys、Zs改正后的值确定初始值c=w=k=0，Xs0，Ys0，Zs0小于跳出循环，完成迭代计算，精度评定判断改正数是否小于限差?大于输入GCP的像点坐标x，y正数小于0.01m,角度改正数小于0.0003，相当于1’的角度值）为止。在这个过程中采用迭代计算的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。2、前方交会七、实验原理公式1、后方交会中运用的共线方程数学模型ZYfZZcYYbXXaZZcYYbXXafyyZXfZZcYYbXXaZZcYYbXXafxxssssssssssss)()()()()()()()()()()()(333222033311103、前方交会与后方交会中均用到旋转矩阵进行的坐标转换输入所需计算点的像平面坐标x,y根据后方交会所得的旋转矩阵Ra,Rb计算像点在左、右像空间辅助坐标系中的坐标XaYaZa，XbYbZb计算摄影基线的三个坐标分量BxByBz计算个点在左右像片中的的投影系数Na,Nb计算地面所求点在地面摄影测量坐标系中的坐标XYZ计算完毕ssssssZZYYXXZZYYXXcbacbacbaZYX1333222111R4、精度评定中均采用最小二乘准则进行平差计算yyZZyYYyXXyyyyvxxZZxYYxXXxxxxvssssssyssssssx00yyaaaZaYaXavxxaaaZaYaXavsssysssx0262524232221016151413121126252423222116151413121100,,aaaaaaaaaaaayyxxZYXvvsssyxAlxVlAxV)()(T1TlAAAx62T0nVV1T)(AAQxxiixximQ05.前方交会的转换2'12'12'1222222111111212121121212,,,ZNBZNZZZYNBYNYYYXNBXNXXXfyxRZYXfyxRZYXZZBYYBXXBZSSYSSXSSSSZSSYSSX标）：地面坐标（摄影测量坐像空间辅助坐标：基线分量：12211121221221ZXZXXBZBNZXZXXBZBNZXZX八．实验源程序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:4X(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);endA=[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')*LV=A*Xp-LXs0=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))breakendend%精度评定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,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);y