1、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计
2、算像点在像空间坐标系中的近似值（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 计算
3、摄影基线的三个坐标分量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
4、（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）; %计算单位权中误
5、差并将其输出显示输入GCP像点坐标及地面摄影测量坐标系坐标的函数和输入所求点像点坐标函数：function shurujcp（）global c m;m=input（摄影比例尺： %输入GCP像点坐标数据函数并分别将其c=input（GCP的总数= % 存入到不同的矩阵之中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 GCP右片像框标坐标：global g2;g2=zeros（c,2）; g2（i,1）=m;g2（i,2）=n; %fun
6、ction shuru（） global a;a=input（计算总像对点数= %完成想计算所需的像平面坐标global b1; %坐标输入，存入不同的矩阵中b1=zeros（a,2）;左片像点坐标：）=a b1（i,1）=m;b1（i,2）=n;end %global b2;b2=zeros（a,2）;右片像点坐标： b2（i,1）=m;b2（i,2）=n;global c;GCP摄影测量系坐标：global G;G=zeros（3,c）;X=Y= v=input（Z= G（i,1）=m;G（i,2）=n;G（i,3）=v;%空间前方交会和后方交会函数：function XYZ=MQZqia
7、nfangjh（）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 %将计算所得的外方位元素存
8、入到.txt % 文件中c g1（i,1）=g1（i,1）+V1（2*i-1）;g1（i,2）=g1（i,2）+V1（2*i）;左片像点坐标.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;右片外方位元素为.txt WF ascii %将计算所得的外方位元素存入到.txt % 文件中c g2（i,1）=g2（i,1）+V2（2*i-1）;g2（i,2）=g2（i,
9、2）+V2（2*i）;右片像点坐标.txt g2 -ascii X1=zeros（a,3）;X2=zeros（a,3）;xx=zeros（3,1）;xxx=zeros（3,1）;a ss=b1（i,1）;b1（i,2）;-f;dd=b2（i,1）;b2（i,2）; xx=R1*ss; X1（i,:）=xx;xxx=R2*dd;X2（i,:）=xxxglobal 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）; N1（1,i）=（BX*X2（i,3
10、）-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）; 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; 左片外方位元素为:Xs Ys Zs disp（WF1
11、）;左片外方位元素协因素阵为：disp（QXX1）;左片像点坐标为：disp（g1）右片外方位元素为:disp（WF2）;右片外方位元素协因素阵为：disp（QXX2）右片像点坐标为：disp（g2）计算所得点摄影测量坐标（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
12、;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
13、_,b3;c1,c2,c3; aX=;aY=;aZ=; for i=1: 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=; xj（i）=-f*aX（i）/aZ（i）; yj（i）=-f*aY（i）/aZ（i）; a11=;a12=;a13=;a14=;a15=;a16=; a21=;a22=;a23=;a2
14、4=;a25=;a26=; 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
15、,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）; lx=;ly=;c lx（i）=g1（i,1）-xj（i）;ly（i）=g1（i,2）-yj（i）; A=zeros（2*c,6）; A（2*i-1,1）=a11（i）;A（2
16、*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）; L=zeros（2*c,1）; L（2*i-1,1）=lx（i）; L（2*i,1）=ly（i）; X=inv（A）*A）*（A）*L; Xs=Xs+X（1,1）;Ys=Ys+X（2,1）;Zs=Zs+X（3,1）;q=
17、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; %后跳出循环，计算结果已经收敛global V;V=Lglobal 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.0001950e+003 5.0007250e+003 2.0201583e+003 -7.2888190e-005 2.819
18、3877e-002 9.5130388e-002 4.0166895e-008 -3.7263703e-010 1.3218695e-008 7.0720033e-005 1.0001730e-007 -2.5748604e-006 -3.7263703e-010 4.0032797e-008 2.6568407e-009 -2.1103715e-007 7.7772275e-005 1.9993587e-005 1.3218695e-008 2.6568407e-009 1.7931301e-008 3.1008915e-005 6.6697659e-006 5.6403374e-007 7
19、.0720033e-005 -2.1103715e-007 3.1008915e-005 1.3087511e-001 1.0148977e-003 -1.9981396e-003 1.0001730e-007 7.7772275e-005 6.6697659e-006 1.0148977e-003 1.5539404e-001 3.0264331e-002 -2.5748604e-006 1.9993587e-005 5.6403374e-007 -1.9981396e-003 3.0264331e-002 4.0721943e-0025.8967023e+003 5.0687355e+00
20、3 2.0506347e+003 1.4337709e-002 4.6257617e-002 1.1037952e-001 3.9305329e-008 4.9400147e-010 -1.0339207e-008 6.8065940e-005 -4.2504770e-007 1.8461496e-006 4.9400147e-010 3.9051893e-008 3.3958896e-011 -3.9945442e-008 7.6312421e-005 -1.6453951e-005 -1.0339207e-008 3.3958896e-011 1.5155886e-008 -2.37050
21、97e-005 3.5940467e-007 -7.3527082e-007 6.8065940e-005 -3.9945442e-008 -2.3705097e-005 1.2229164e-001 -2.3449223e-003 4.8281474e-003 -4.2504770e-007 7.6312421e-005 3.5940467e-007 -2.3449223e-003 1.5233230e-001 -2.5374659e-002 1.8461496e-006 -1.6453951e-005 -7.3527082e-007 4.8281474e-003 -2.5374659e-0
22、02 3.6794789e-002GCP在左片和右片改正后的坐标（x，y）为： 1.6019582e-002 7.9954660e-002 -7.3934212e-002 7.8699356e-002 8.8559633e-002 8.1141190e-002 -5.2455612e-003 7.8187184e-002 1.3352398e-002 -7.9378247e-002 -7.9125440e-002 -7.8877760e-002 8.2242309e-002 -8.0017749e-002 -9.8858970e-003 -8.0086832e-002单位权中误差为：1.515610577029578e-005 所求地面点的三维坐标（X, Y, Z）为： 5.4310348e+003 5.8851463e+003 5.4831646e+002 5.1473645e+003 5.0555934e+003 4.8499600e+002 5.4957931e+003 5.0826911e+003 5.0668967e+002 5.8442434e+003 5.1098033e+003 5.3025650e+002 5.5603279e+003 4.2870779e+003 4.6536459e+0
