摄影测量实验报告(空间后方交会前方交会)Word格式.docx
《摄影测量实验报告(空间后方交会前方交会)Word格式.docx》由会员分享，可在线阅读，更多相关《摄影测量实验报告(空间后方交会前方交会)Word格式.docx（11页珍藏版）》请在冰点文库上搜索。
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的总数='
%存入到不同的矩阵之中
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
GCP右片像框标坐标：
globalg2;
g2=zeros（c,2）;
g2（i,1）=m;
g2（i,2）=n;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionshuru（）
globala;
a=input（'
计算总像对点数='
%完成想计算所需的像平面坐标
globalb1;
%坐标输入，存入不同的矩阵中
b1=zeros（a,2）;
左片像点坐标：
）
=a
b1（i,1）=m;
b1（i,2）=n;
end
%%
globalb2;
b2=zeros（a,2）;
右片像点坐标：
b2（i,1）=m;
b2（i,2）=n;
globalc;
GCP摄影测量系坐标：
globalG;
G=zeros（3,c）;
X='
Y='
v=input（'
Z='
G（i,1）=m;
G（i,2）=n;
G（i,3）=v;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
空间前方交会和后方交会函数：
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
%文件中
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,k2R]=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,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,:
）=xxx'
globalXs1Xs2Ys1Ys2Zs1Zs2;
BX=Xs2-Xs1;
BY=Ys2-Ys1;
BZ=Zs2-Zs1;
globalN1N2;
N1=zeros（1,a）;
N2=zeros（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）;
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;
左片外方位元素为:
XsYsZsψωκ'
disp（WF1）;
左片外方位元素协因素阵为：
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,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:
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=[];
a24=[];
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,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*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=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;
%后跳出循环，计算结果已经收敛
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
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+0
升级会员 