摄影测量实验报告材料空间后方交会前方交会.docx
《摄影测量实验报告材料空间后方交会前方交会.docx》由会员分享，可在线阅读，更多相关《摄影测量实验报告材料空间后方交会前方交会.docx（18页珍藏版）》请在冰豆网上搜索。
摄影测量实验报告材料空间后方交会前方交会
空间后方交会-空间前方交会程序编程实验
一．实验目的要求
掌握运用空间后方交会-空间前方交会求解地面点的空间位置。
学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。
然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程。
二．仪器用具
计算机、编程软件（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
组成法方程式，计算系数A’A常数项A’L求解外方位元素
计算ψ、ω、κ、Xs、Ys、Zs改正后的值
大于
小于
计算完毕
此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0.00003，相当于0.1’的角度值）为止。
在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。
输入所需计算点的像平面坐标x1,y1；x2,y2
根据后方交会所得的旋转矩阵R1,R2计算像点在左右像空间辅助坐标系中的坐标X1Y1Z1，X2Y2Z2
计算摄影基线的三个坐标分量BxByBz
计算个点在左右像片中的的投影系数N1N2
计算地面所求点在地面摄影测量坐标系中的坐标
XAYAZA
计算完毕，精度评定
在空间后方交会中运用的数学模型为共线方程
确定Xs，Ys，Zs的初始值时，对于左片可取地面左边两个GCP的坐标的平均值作为左片Xs和Ys的初始值，取右边两个GCP的坐标平均值作为右片Xs和Ys的初始值。
Zs可取地面所有GCP的Z坐标的平均值再加上航高。
空间前方交会的数学模型为：
五．实验源代码
functionMain_KJQHFJH（）
globalRg1g2mGacb1b2;
m=10000;a=5;c=4;
feval（@shuru）;%调用shuru（）shurujcp（）函数完成像点及
feval（@shurujcp）;%CCP有关数据的输入
XYZ=feval（@MQZqianfangjh）;%调用MQZqianfangjh（）函数完成空间前方、
%%%%%%单位权中误差%%%%%后方交会计算解得外方位元素
globalV1V2;%由于以上三个函数定义在外部文件中故需
VV=[];%用feval（）完成调用过程
fori=1:
2*c
VV（i）=V1（i）;VV（2*i+1）=V2（i）;
end
m0=sqrt（VV*（VV'）/（2*c-6））;
disp（'单位权中误差m0为正负:
'）;disp（m0）;%计算单位权中误差并将其输出显示
输入GCP像点坐标及地面摄影测量坐标系坐标的函数和输入所求点像点坐标函数：
functionshurujcp（）
globalcm;
m=input（'摄影比例尺：
'）;%输入GCP像点坐标数据函数并分别将其
c=input（'GCP的总数='）;%存入到不同的矩阵之中
disp（'GCP左片像框标坐标：
'）;
globalg1;g1=zeros（c,2）;
i=1;
whilei<=c
m=input（'x='）;
n=input（'y='）;
g1（i,1）=m;g1（i,2）=n;
i=i+1;
end
disp（'GCP右片像框标坐标：
'）;
globalg2;g2=zeros（c,2）;
i=1;
whilei<=c
m=input（'x='）;
n=input（'y='）;
g2（i,1）=m;g2（i,2）=n;
i=i+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionshuru（）
globala;
a=input（'计算总像对点数='）;%完成想计算所需的像平面坐标
globalb1;%坐标输入，存入不同的矩阵中
b1=zeros（a,2）;
disp（'左片像点坐标：
'）
i=1;
whilei<=a
m=input（'x='）;
n=input（'y='）;
b1（i,1）=m;b1（i,2）=n;
i=i+1;
end
%%
globalb2;
b2=zeros（a,2）;
disp（'右片像点坐标：
'）
i=1;
whilei<=a
m=input（'x='）;
n=input（'y='）;
b2（i,1）=m;b2（i,2）=n;
i=i+1;
end
%%
globalc;
c=input（'GCP的总数='）;
disp（'GCP摄影测量系坐标：
'）
globalG;
G=zeros（3,c）;
i=1;
whilei<=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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
空间前方交会和后方交会函数：
functionXYZ=MQZqianfangjh（）
globalR1R2afb1b2RaRb;
globalX1X2;R1=Ra;R2=Rb;R1=zeros（3,3）;R2=zeros（3,3）;
globalg1g2GV1V2VWFcQXXQXX1QXX2;
xs0=（G（1,1）+G（3,1））/2;
ys0=（G（1,2）+G（3,2））/2;
[Xs1,Ys1,Zs1,q1,w1,k1R]=houfangjh（g1,xs0,ys0）;%对左片调用后方交会函数
R1=R;
V1=V;
WF1=WF;
QXX1=QXX;
save'左片外方位元素为.txt'WF-ascii%将计算所得的外方位元素存入到.txt
%文件中
fori=1:
c
g1（i,1）=g1（i,1）+V1（2*i-1）;g1（i,2）=g1（i,2）+V1（2*i）;
end
save'左片像点坐标.txt'g1-ascii
xs0=（G（2,1）+G（4,1））/2;
ys0=（G（2,2）+G（4,2））/2;
[Xs2,Ys2,Zs2,q2,w2,k2R]=houfangjh（g2,xs0,ys0）;%对右片调用后方交会函数
R2=R;
V2=V;
WF2=WF;
QXX2=QXX;
save'右片外方位元素为.txt'WF–ascii%将计算所得的外方位元素存入到.txt
%文件中
fori=1:
c
g2（i,1）=g2（i,1）+V2（2*i-1）;g2（i,2）=g2（i,2）+V2（2*i）;
end
save'右片像点坐标.txt'g2-ascii
X1=zeros（a,3）;X2=zeros（a,3）;
xx=zeros（3,1）;xxx=zeros（3,1）;
fori=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';
end
globalXs1Xs2Ys1Ys2Zs1Zs2;
BX=Xs2-Xs1;BY=Ys2-Ys1;BZ=Zs2-Zs1;
globalN1N2;N1=zeros（1,a）;N2=zeros（1,a）;
fori=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%计算投影系数，并计算五点的三维坐标
globalXYZ;XYZ=zeros（a,3）;
fori=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;
end
disp（'左片外方位元素为:
XsYsZsψωκ'）;
disp（WF1）;
disp（'左片外方位元素协因素阵为：
'）;
disp（QXX1）;
disp（'左片像点坐标为：
'）
disp（g1）
disp（'右片外方位元素为:
XsYsZsψωκ'）;
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,kR]=houfangjh（g1,Xs0,Ys0）%计算像片外方位元素
%%%%%%%%%%
globalfGmcb1b2;
f=0.152;
Xs=Xs0;
Ys=Ys0;
Zs=m*f+G（1,3）;
q=0;w=0;k=0;
while1%实现一个永真循环，是改正数小于限差以后跳出循环
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=[];
fori=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=[];
fori=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=[];
fori=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=[];
fori=1:
c
lx（i）=g1（i,1）-xj（i）;ly（i）=g1（i,2）-yj（i）;
end
A=zeros（2*c,6）;
fori=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）;
fori=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）;
ifaaa<0.00003%当改正数中绝对值最大的改正数小于限差0.00003
break;%后跳出循环，计算结果已经收敛
end
end
globalV;V=L';
globalWFQXX;
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+0035.0007250e+0032.0201583e+003-7.2888190e-0052.8193877e-0029.5130388e-002
左片外方位元素协因素阵为：
4.0166895e-008-3.7263703e-0101.3218695e-0087.0720033e-0051.0001730e-007-2.5748604e-006
-3.7263703e-0104.0032797e-0082.6568407e-009-2.1103715e-0077.7772275e-0051.9993587e-005
1.3218695e-0082.6568407e-0091.7931301e-0083.1008915e-0056.6697659e-0065.6403374e-007
7.0720033e-005-2.1103715e-0073.1008915e-0051.3087511e-0011.0148977e-003-1.9981396e-003
1.0001730e-0077.7772275e-0056.6697659e-0061.0148977e-0031.5539404e-0013.0264331e-002
-2.5748604e-0061.9993587e-0055.6403374e-007-1.9981396e-0033.0264331e-0024.0721943e-002
左片外方位元素Xs,Ys,Zs,ψ、ω、κ、为：
5.8967023e+0035.0687355e+0032.0506347e+0031.4337709e-0024.6257617e-0021.1037952e-001
右片外方位元素协因素阵为：
3.9305329e-0084.9400147e-010-1.0339207e-0086.8065940e-005-4.2504770e-0071.8461496e-006
4.9400147e-0103.9051893e-0083.3958896e-011-3.9945442e-0087.6312421e-005-1.6453951e-005
-1.0339207e-0083.3958896e-0111.5155886e-008-2.3705097e-0053.5940467e-007-7.3527082e-007
6.8065940e-005-3.9945442e-008-2.3705097e-0051.2229164e-001-2.3449223e-0034.8281474e-003
-4.2504770e-0077.6312421e-0053.5940467e-007-2.3449223e-0031.5233230e-001-2.5374659e-002
1.8461496e-006-1.6453951e-005-7.3527082e-0074.8281474e-003-2.5374659e-0023.6794789e-002
GCP在左片和右片改正后的坐标（x，y）为：
1.6019582e-0027.9954660e-002-7.3934212e-0027.8699356e-002
8.8559633e-0028.1141190e-002-5.2455612e-0037.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+0035.8851463e+0035.4831646e+002
5.1473645e+0035.0555934e+0034.8499600e+002
5.4957931e+0035.0826911e+0035.0668967e+002
5.8442434e+0035.1098033e+0035.3025650e+002
5.5603279e+0034.2870779e+0034.6536459e+002
七．心得体会
经过三周的努力，这个当初看来艰巨的任务终于在我的不懈努力下圆满的完成了。
在编程的过程中，不仅加强了我对空间后方交会求解外方位元素以及应用空间前方交会求解物点在摄影测量坐标系中的三维坐标，同时也锻炼了我们将理论知识应用到实际中的一种能力。
在这个过程中我遇到了很多的困难，经过查看有关的资料，并在同学的帮助下成功的得到了解决。
在这之前我们并没有系统的学过MATLAB这个强大的数学处理软件，对它了解的太少，但我还是毅然的选择了它。
因为对于MATLAB，它在处理有关矩阵运算的时候显示出了太强大的优势，与其他编程语言来说，简单易操作。
在老师刚刚把任务交给我们的时候，我并没与急切的就去实践。
而是将空间后方交会和空间前方交会有关的理论知识系统的看了几遍。
在上课的时候，老师一讲好像什么都明白了似的，但是每当合上课本的时候，好像什么都忘记了，为了能有一个比较完整的思路，我一边看书上介绍的求解方法，一边就把具体的求解过程用流程图在草稿上粗略的表达了一下，这对我后面的编程过程有极大的帮助作用。
在对书本知识有了一个系统完善的理解以后，接下了就是一个艰巨的过程了，对于理论知识我们我懂我们可以多看几遍就可以有一个粗略的理解，但对于实践来说，首先我们必须对理论有一个透彻的理解，然后还要将它转化的实际应用中去。
在刚刚开始编写代码的时候，只是慢慢的不断地摸索着，对于一个语句的定义，一个结构的格式我们都只能到书本和网络上一步一步的求解。
这在刚刚开始的时候可以说对我来说是一个极大的考验。
经过了几天的临时突击，对于编程的语言我已经有了一个大体的了解，在整个程序中所需设计到得，我基本都已经掌握了。
接下来就是，具体功能的编写了。
根据前面所总结的流程图，我对整个程序的框架有了一个比较完善的把握了。
整个大的程序中，我将它们分割成为几个小的功能，对于每个功能，一步一步的来实现。
整个程序我将它分成了三大块，第一就是数据的输入，第二就是空间后方交会求解像片的外方位元素，第三就是空间前方交会求解地面点的三维坐标。
在完成了有关代码的编写以后，就是程序的调试了，这个看似比较容易的过程，然而却花了我最多的时间。
在刚一开始一运行，全是错误，我就只能一句一句
升级会员 