摄影测量实验报告(空间后方交会—前方交会)

空间后方交会-空间前方交会程序编程实验一 实验目的要求掌握运用空间后方交会-空间前方交会求解地面点的空间位置。学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程。二 仪器用具 计算机、编程软件（MATLAB）三 实验数据实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。此四对坐标运用最小二乘法求解左右像片的外方位元素，即完成了空间后方的过程。另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。实验数据如下：点号 左片 右片 地面摄影测量坐标 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 81.555 -39.953 78.463 2 14.618 -0.231 -76.016 0.036 3 49.88 -0.792 -42.201 -1.022 4 86.243 -1.346 -7.706 -2.112 5 48.135 -79.962 -44.438 -79.736 内方位元素：f=152.000mm，x0=0，y0=0四 实验框图输入GCP的像点坐标 xy确定初始值=0，Xs，Ys，Zs计算旋转矩阵R计算像点在像空间坐标系中的近似值（x）,（y）,并组成误差方程的常数项l计算误差方程的系数项组成系数矩阵A组成法方程式，计算系数AA常数项AL求解外方位元素计算、Xs、Ys、Zs改正后的值改正数是否小于限差？ 大于小于计算完毕此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0.00003，相当于0.1的角度值）为止。在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。输入所需计算点的像平面坐标x1,y1；x2,y2根据后方交会所得的旋转矩阵R1,R2计算像点在左右像空间辅助坐标系中的坐标X1 Y1 Z1，X2 Y2 Z2 计算摄影基线的三个坐标分量Bx By Bz 计算个点在左右像片中的的投影系数N1 N2计算地面所求点在地面摄影测量坐标系中的坐标XA YA ZA计算完毕，精度评定在空间后方交会中运用的数学模型为共线方程确定Xs，Ys，Zs的初始值时，对于左片可取地面左边两个GCP的坐标的平均值作为左片Xs 和Ys的初始值，取右边两个GCP的坐标平均值作为右片Xs 和Ys的初始值。Zs可取地面所有GCP的Z坐标的平均值再加上航高。空间前方交会的数学模型为： 五 实验源代码function Main_KJQHFJH()global R g1 g2 m G a c b1 b2;m=10000;a=5;c=4;feval(shuru); %调用shuru（）shurujcp（）函数完成像点及feval(shurujcp); %CCP有关数据的输入 XYZ=feval(MQZqianfangjh); %调用MQZqianfangjh（）函数完成空间前方、% 单位权中误差 % %后方交会计算解得外方位元素global V1 V2; %由于以上三个函数定义在外部文件中故需VV=; %用feval（）完成调用过程for i=1:2*c VV(i)=V1(i);VV(2*i+1)=V2(i);endm0=sqrt(VV*(VV)/(2*c-6);disp(单位权中误差m0为正负:);disp(m0); %计算单位权中误差并将其输出显示输入GCP像点坐标及地面摄影测量坐标系坐标的函数和输入所求点像点坐标函数：function shurujcp()global c m;m=input(摄影比例尺：); %输入GCP像点坐标数据函数并分别将其c=input(GCP的总数=); % 存入到不同的矩阵之中disp(GCP左片像框标坐标：);global g1;g1=zeros(c,2);i=1;while i=c m=input(x=); n=input(y=); g1(i,1)=m;g1(i,2)=n; i=i+1;end disp(GCP右片像框标坐标：);global g2;g2=zeros(c,2);i=1;while i=c m=input(x=); n=input(y=); g2(i,1)=m;g2(i,2)=n; i=i+1; end%function shuru() global a;a=input(计算总像对点数=); %完成想计算所需的像平面坐标global b1; %坐标输入，存入不同的矩阵中b1=zeros(a,2);disp(左片像点坐标：)i=1;while i=a m=input(x=); n=input(y=); b1(i,1)=m;b1(i,2)=n; i=i+1;end %global b2;b2=zeros(a,2);disp(右片像点坐标：)i=1;while i=a m=input(x=); n=input(y=); b2(i,1)=m;b2(i,2)=n; i=i+1;end %global c;c=input(GCP的总数=);disp(GCP摄影测量系坐标：)global G;G=zeros(3,c);i=1;while i=c m=input(X=); n=input(Y=); v=input(Z=); G(i,1)=m;G(i,2)=n;G(i,3)=v; i=i+1;end%空间前方交会和后方交会函数：function XYZ=MQZqianfangjh()global R1 R2 a f b1 b2 Ra Rb;global X1 X2;R1=Ra;R2=Rb;R1=zeros(3,3);R2=zeros(3,3);global g1 g2 G V1 V2 V WF c QXX QXX1 QXX2;xs0=(G(1,1)+G(3,1)/2;ys0=(G(1,2)+G(3,2)/2;Xs1,Ys1,Zs1,q1,w1,k1 R=houfangjh(g1,xs0,ys0); %对左片调用后方交会函数R1=R;V1=V;WF1=WF;QXX1=QXX;save 左片外方位元素为.txt WF -ascii %将计算所得的外方位元素存入到.txt % 文件中for i=1:c g1(i,1)=g1(i,1)+V1(2*i-1);g1(i,2)=g1(i,2)+V1(2*i);endsave 左片像点坐标.txt g1 -ascii xs0=(G(2,1)+G(4,1)/2;ys0=(G(2,2)+G(4,2)/2;Xs2,Ys2,Zs2,q2,w2,k2 R=houfangjh(g2,xs0,ys0); %对右片调用后方交会函数 R2=R; V2=V;WF2=WF;QXX2=QXX;save 右片外方位元素为.txt WF ascii %将计算所得的外方位元素存入到.txt % 文件中for i=1:c g2(i,1)=g2(i,1)+V2(2*i-1);g2(i,2)=g2(i,2)+V2(2*i);endsave 右片像点坐标.txt g2 -ascii X1=zeros(a,3);X2=zeros(a,3);xx=zeros(3,1);xxx=zeros(3,1);for i=1:a ss=b1(i,1);b1(i,2);-f;dd=b2(i,1);b2(i,2);-f; xx=R1*ss; X1(i,:)=xx;xxx=R2*dd;X2(i,:)=xxx;endglobal Xs1 Xs2 Ys1 Ys2 Zs1 Zs2; BX=Xs2-Xs1;BY=Ys2-Ys1;BZ=Zs2-Zs1;global N1 N2;N1=zeros(1,a);N2=zeros(1,a);for i=1:a N1(1,i)=(BX*X2(i,3)-BZ*X2(i,1)/(X1(i,1)*X2(i,3)-X2(i,1)*X1(i,3); N2(1,i)=(BX*X1(i,3)-BZ*X1(i,1)/(X1(i,1)*X2(i,3)-X2(i,1)*X1(i,3);end %计算投影系数，并计算五点的三维坐标global XYZ;XYZ=zeros(a,3);for i=1:a XYZ(i,1)=Xs1+N1(1,i)*X1(i,1); XYZ(i,3)=Zs1+N1(1,i)*X1(i,3); XYZ(i,2)=(Ys1+N1(1,i)*X1(i,2)+(Ys2+N2(1,i)*X2(i,2)/2; enddisp(左片外方位元素为:Xs Ys Zs );disp(WF1);disp(左片外方位元素协因素阵为：);disp(QXX1);disp(左片像点坐标为：)disp(g1)disp(右片外方位元素为:Xs Ys Zs );disp(WF2);disp(右片外方位元素协因素阵为：)disp(QXX2)disp(右片像点坐标为：)disp(g2)disp(计算所得点摄影测量坐标（X,Y,Z）为：);disp(XYZ); save XYZ.txt XYZ -ascii %将计算所得结果保存到XYZ.txt文件中%function Xs,Ys,Zs,q,w,k R=houfangjh(g1,Xs0,Ys0) %计算像片外方位元素%global f G m c b1 b2;f=0.152; Xs=Xs0;Ys=Ys0;Zs=m*f+G(1,3);q=0;w=0;k=0;while 1 %实现一个永真循环，是改正数小于限差以后跳出循环 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; aX=;aY=;aZ=; for i=1:c aX(i)=a1*(G(i,1)-Xs)+b1_*(G(i,2)-Ys)+c1*(G(i,3)-Zs); aY(i)=a2*(G(i,1)-Xs)+b2_*(G(i,2)-Ys)+c2*(G(i,3)-Zs); aZ(i)=a3*(G(i,1)-Xs)+b3*(G(i,2)-Ys)+c3*(G(i,3)-Zs); end xj=;yj=; for i=1:c xj(i)=-f*aX(i)/aZ(i); yj(i)=-f*aY(i)/aZ(i); end a11=;a12=;a13=;a14=;a15=;a16=; a21=;a22=;a23=;a24=;a25=;a26=; for i=1:c a11(i)=(a1*f+a3*g1(i,1)/aZ(i);a12(i)=(b1_*f+b3*g1(i,1)/aZ(i);a13(i)=(c1*f+c3*g1(i,1)/aZ(i); a21(i)=(a2*f+a3*g1(i,2)/aZ(i);a22(i)=(b2_*f+b3*g1(i,2)/aZ(i);a23(i)=(c2*f+c3*g1(i,2)/aZ(i); a14(i)=g1(i,2)*sin(w)-(g1(i,1)*(g1(i,1)*cos(k)-g1(i,2)*sin(k)/f+f*cos(k)*cos(w); a15(i)=-f*sin(k)-g1(i,1)*(g1(i,1)*sin(k)+g1(i,2)*cos(k)/f; a16(i)=g1(i,2); a24(i)=-g1(i,1)*sin(w)-(g1(i,2)*(g1(i,1)*cos(k)-g1(i,2)*sin(k)/f-f*sin(k)*cos(w); a25(i)=-f*cos(k)-g1(i,2)*(g1(i,1)*sin(k)+g1(i,2)*cos(k)/f; a26(i)=-g1(i,1); end lx=;ly=;for i=1:c lx(i)=g1(i,1)-xj(i);ly(i)=g1(i,2)-yj(i); end A=zeros(2*c,6); for i=1:c A(2*i-1,1)=a11(i);A(2*i-1,2)=a12(i);A(2*i-1,3)=a13(i);A(2*i-1,4)=a14(i);A(2*i-1,5)=a15(i);A(2*i-1,6)=a16(i); A(2*i,1)=a21(i); A(2*i,2)=a22(i); A(2*i,3)=a23(i); A(2*i,4)=a24(i); A(2*i,5)=a25(i); A(2*i,6)=a26(i); end L=zeros(2*c,1); for i=1:c L(2*i-1,1)=lx(i); L(2*i,1)=ly(i); end X=inv(A)*A)*(A)*L; Xs=Xs+X(1,1);Ys=Ys+X(2,1);Zs=Zs+X(3,1);q=q+X(4,1);w=w+X(5,1);k=k+X(6,1); Xabs=abs(X); aaa=max(Xabs); if aaa0.00003 %当改正数中绝对值最大的改正数小于限差0.00003 break; %后跳出循环，计算结果已经收敛 endendglobal V;V=L;global WF QXX;WF(1)=Xs;WF(2)=Ys;WF(3)=Zs;WF(4)=q;WF(5)=w;WF(6)=k;QXX=A*A;六 实验结果左片外方位元素Xs,Ys,Zs, 、为：5.e+003 5.e+003 2.e+003 -7.e-005 2.e-002 9.e-002左片外方位元素协因素阵为： 4.e-008 -3.e-010 1.e-008 7.e-005 1.e-007 -2.e-006 -3.e-010 4.e-008 2.e-009 -2.e-007 7.e-005 1.e-005 1.e-008 2.e-009 1.e-008 3.e-005 6.e-006 5.e-007 7.e-005 -2.e-007 3.e-005 1.e-001 1.e-003 -1.e-003 1.e-007 7.e-005 6.e-006 1.e-003 1.e-001 3.e-002 -2.e-006 1.e-005 5.e-007 -1.e-003 3.e-002 4.e-002左片外方位元素Xs,Ys,Zs, 、为： 5.e+003 5.e+003 2.e+003 1.e-002 4.e-002 1.e-001右片外方位元素协因素阵为： 3.e-008 4.e-010 -1.e-008 6.e-005 -4.e-007 1.e-006 4.e-010 3.e-008 3.e-011 -3.e-008 7.e-005 -1.e-005 -1.e-008 3.e-011 1.e-008 -2.e-005 3.e-007 -7.e-007 6.e-005 -3.e-008 -2.e-005 1.e-001 -2.e-003 4.e-003 -4.e-007 7.e-005 3.e-007 -2.e-003 1.e-001 -2.e-002 1.e-006 -1.e-005 -7.e-007 4.e-003 -2.e-002 3.e-002GCP在左片和右片改正后的坐标（x，y）为： 1.e-002 7.e-002 -7.e-002 7.e-002 8.e-002 8.e-002 -5.e-003 7.e-002 1.e-002 -7.e-002 -7.e-002 -7.e-002 8.e-002 -8.e-002 -9.e-003 -8.e-002单位权中误差为：1.9578e-005 所求地面点的三维坐标（X, Y, Z）为： 5.e+003 5.e+003 5.e+002 5.e+003 5.e+003 4.e+002 5.e+003 5.e+003 5.e+002 5.e+003 5.e+003 5.e+002 5.e+003 4.e+003 4.e+002七 心得体会经过三周的努力，这个当初看来艰巨的任务终于在我的不懈努力下圆满的完成了。在编程的过程中，不仅加强了我对空间后方交会求解外方位元素以及应用空间前方交会求解物点在摄影测量坐标系中的三维坐标，同时也锻炼了我们将理论知识应用到实际中的一种能力。在这个过程中我遇到了很多的困难，经过查看有关的资料，并在同学的帮助下成功的得到了解决。在这之前我们并没有系统的学过MATLAB这个 强大的数学处理软件，对它了解的太少，但我还是毅然的选择了它。因为对于MATLAB，它在处理有关矩阵运算的时候显示出了太强大的优势，与其他编程语言来说，简单易操作。在老师刚刚把任务交给我们的时候，我并没与急切的就去实践。而是将空间后方交会和空间前方交会有关的理论知识系统的看了几遍。在上课的时候，老师一讲好像什么都明白了似的，但是每当合上课本的时候，好像什么都忘记了，为了能有一个比较完整的思路，我一边看书上介绍的求解方法，一边就把具体的求解过程用流程图在草稿上粗略的表达了一下，这对我后面的编程过程有极大的帮助作用。在对书本知识有了一个系统完善的理解以后，接下了就是一个艰巨的过程了，对于理论知识我们我懂我们可以多看几遍就可以有一个粗略的理解，但对于实践来说，首先我们必须对理论有一个透彻的理解，然后还要将它转化的实际应用中去。在刚刚开始编写代码的时候，只是慢慢的不断地摸索着，对于一个语句的定义，一个结构的格式我们都只能到书本和网络上一步一步的求解。这在刚刚开始的时