摄影测量学
实验报告
学院： 地信院
班级： 测绘0904班
老师： 邹峥嵘
姓名： 张文佳
学号： 0405090921
20##年11 月11 日
空间后方交会——空间前方交会
程序编程实验
一． 实验目的
1、要求掌握运用摄影测量中空间后方交会-空间前方交会求解地面点的空间位置的方法和原理。
2、学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的计算，完成所给像对中两张像片各自的六个外方位元素的求解和精度评定。
3、根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，利用计算机编程语言前方交会编程，求解其对应的地面点在摄影测量坐标系中的坐标，从而达到通过摄影测量量测地面地理数据的目的。
二． 实验仪器
1、计算机
2、MATLAB计算机编程软件
三、实验数据
实验数据实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。此四对坐标运用最小二乘法求解左右像片的外方位元素， 即完成了空间后方的过程。另外还 给出了 5 对地面点在左右像片中的像平面坐标和左右像片的内方位元素。实验数据如下：
f=150.000mm，x0=0，y0=0
四、程序设计流程图
1、后方交会
此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（长度改正数小于0.01m,角度改正数小于0.0003，相当于1’的角度值）为止。在这个过程中采用迭代计算的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。
2、前方交会
七、实验原理公式
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 %统计迭代计算次数
%完成迭代计算，检验改正是是否符合要求
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
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
for 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)));
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.0951 0.9950 -0.0291
0.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+003
Ys0 = 5.0001e+003
Zs0 = 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+003
Ys0 = 5.0879e+003
Zs0 = 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.5274
Y =
1.0e+003 *
5.7446 5.0165 5.0414 5.0679 4.2922
Z =
1.0e+003 *
1.8901 1.8638 1.8633 1.8597 1.8368
十、实验心得体会
经过几个星期的努力，这个当初看似不可能完成的任务终于完成了，我感到很欣慰很有成就感。通过这次编程实验，加强了我对空间后方交会求解外方元素以及应用空间前方交会求解物点在摄影测量做坐标系中的三维坐标，同时让我认识到将理论联系实际的重要性。
在本次试验中，我都遇到了很多困难，但都被我逐级突破。在这之前我并没有系统性的学习过MATLAB,甚至可以说是根本没用过，虽然通过了计算机二级，但如果仅凭自己的能力，用C编写这个程序，还是觉得毫无头绪。后来听同学说MATLAB简单易学，容易掌握，而且处理矩阵运算有强大的优势，所以我毅然选择了MATLAB作为本次实验编程的软件。花了几天的时间学习了MATLAB的编绘语言以及运算法则。
刚开始编写时毫无头绪，而且觉得变量多得错综复杂，搞得我两眼发晕。后来我没有急着去编写，而是回归到课本，将空间后方交会和前方交会有关的理论知识系统认真的看了几遍，了解了整个计算过程后，弄清变量之间的转换关系后，再将整个计算过程绘制成流程图，再进行编程。
完成了有关代码的基本编写后，就开始了程序的调试了，这是让我最为痛苦的一个环节，看似简单，而我却花了很长时间，语法错误比较容易发现，根据软件提示是很容易找出来的，其中很容易出现的错误是错误的拼写，矩阵维数不一致等问题。就在没有语法错误的时候我满怀希望的再次运行，却得到的了无限循环，让我很吃惊，强制退出程序后，重新检查了很多遍源程序，把书上的公式和自己的进行了详细比对，最后才发现因为马虎导致有的公式漏打括号，变量输错等原因，进行改正后得到了正确结果。
在整个过程中，因为时间等原因，我没有对MATLAB进行深入学习，没有运用到程序函数功能，用了两个m文件进行了左片右片的程序编写，导致了很多重复步骤，效率低下，也没有注意到程序运行结果的美观和清晰明了程度。在下次进行运用时会多加注意。
总之，整个实验算是完成了，我不仅对后方交会和前方交会的实用原理更加明了，初步学习了MATLAB,通过本次实验也增强了我分析问题和处理问题的能力，也给了我一个机会把理论知识运用到实际计算中去的机会。