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

空间后方交会-空间前方交会程序编程实验一．实验目的要求

掌握运用空间后方交会-空间前方交会求解地面点的空间位置。学会运用空间后方交会的原理，根据所给控制点的地面摄影测量坐标系坐标以及相应的像平面坐标系中的坐标，利用计算机编程语言实现空间后方交会的过程，完成所给像对中两张像片各自的外方位元素的求解。然后根据空间后方交会所得的两张像片的内外方位元素，利用同名像点在左右像片上的坐标，求解其对应的地面点在摄影测量坐标系中的坐标，并完成精度评定过程，利用计算机编程语言实现此过程。

二．仪器用具

计算机、编程软件（MATLAB）

三．实验数据

实验数据包含四个地面控制点（GCP）的地面摄影测量坐标及在左右像片中的像平面坐标。此四对坐标运用最小二乘法求解左右像片的外方位元素，即完成了空间后方的过程。另外还给出了5对地面点在左右像片中的像平面坐标和左右像片的内方位元素。实验数据如下：

内方位元素：f=152.000mm，x0=0，y0=0

四．实验框图

此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0.00003，相当于0.1’的角度值）为止。在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。

在空间后方交会中运用的数学模型为共线方程

确定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);

end

m0=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);

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,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);

end

save '右片像点坐标.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';

end

global 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;

end

disp('左片外方位元素为: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 aaa<0.00003 %当改正数中绝对值最大的改正数小于限差0.00003 break; %后跳出循环，计算结果已经收敛

end

end

global 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.0001950e+003 5.0007250e+003 2.0201583e+003 -7.2888190e-005 2.8193877e-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.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-002

左片外方位元素Xs,Ys,Zs,ψ、ω、κ、为：

5.8967023e+003 5.0687355e+003 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.3705097e-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-002

3.6794789e-002

GCP在左片和右片改正后的坐标（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+002

七．心得体会

经过三周的努力，这个当初看来艰巨的任务终于在我的不懈努力下圆满的完成了。在编程的过程中，不仅加强了我对空间后方交会求解外方位元素以及应用空间前方交会求解物点在摄影测量坐标系中的三维坐标，同时也锻炼了我们将理论知识应用到实际中的一种能力。在这个过程中我遇到了很多的困难，经过查看有关的资料，并在同学的帮助下成功的得到了解决。

在这之前我们并没有系统的学过MATLAB这个强大的数学处理软件，对它了解的太少，但我还是毅然的选择了它。因为对于MATLAB，它在处理有关矩阵运算的时候显示出了太强大的优势，与其他编程语言来说，简单易操作。在老师刚刚把任务交给我们的时候，我并没与急切的就去实践。而是将空间后方交会和空间前方交会有关的理论知识系统的看了几遍。在

上课的时候，老师一讲好像什么都明白了似的，但是每当合上课本的时候，好像什么都忘记了，为了能有一个比较完整的思路，我一边看书上介绍的求解方法，一边就把具体的求解过程用流程图在草稿上粗略的表达了一下，这对我后面的编程过程有极大的帮助作用。在对书本知识有了一个系统完善的理解以后，接下了就是一个艰巨的过程了，对于理论知识我们我懂我们可以多看几遍就可以有一个粗略的理解，但对于实践来说，首先我们必须对理论有一个透彻的理解，然后还要将它转化的实际应用中去。在刚刚开始编写代码的时候，只是慢慢的不断地摸索着，对于一个语句的定义，一个结构的格式我们都只能到书本和网络上一步一步的求解。这在刚刚开始的时候可以说对我来说是一个极大的考验。经过了几天的临时突击，对于编程的语言我已经有了一个大体的了解，在整个程序中所需设计到得，我基本都已经掌握了。接下来就是，具体功能的编写了。根据前面所总结的流程图，我对整个程序的框架有了一个比较完善的把握了。整个大的程序中，我将它们分割成为几个小的功能，对于每个功能，一步一步的来实现。整个程序我将它分成了三大块，第一就是数据的输入，第二就是空间后方交会求解像片的外方位元素，第三就是空间前方交会求解地面点的三维坐标。

在完成了有关代码的编写以后，就是程序的调试了，这个看似比较容易的过程，然而却花了我最多的时间。在刚一开始一运行，全是错误，我就只能一句一句的来对它进行修改了，在这个过程中，主要是由于粗心而造成的一些语法错误，很快程序基本已经能‘顺利’的往下执行了，但是始终出不来结果，主要是改正数不能收敛，程序已经是一个死循环了。这个问题显然是由于公式的运用错误造成的，在将书本上的数学模型与所编写的公式做了一个比较细致的对比以后，最终成功的找到了问题的所在。

在这次试验中，我只是完成了最基本的要求，由于时间的原因，没有过多的去追求程序的界面，我想在今后的实际中，我们不仅需要编写出高效的程序，而且一个友好的用户界面也是必不可少的。在进行数据交互的时候，在程序的运行过程中输入，容易输错而且但数据量巨大的时候，这种手动输入的方法将难以实现，这也是本程序的一个缺点。对于数据的交互在实际的程序中可以先将其存入到一个txt 获excel 文件中然后在程序的运行过程中直接读取文件中的数据，这样减轻了用户的劳动量，提高了结果的准确信。我想在以后的程序设计方面这也是一个不可忽视的方面。

回顾此处实验，我比较圆满的完成了实验任务，但最关键的还是在这个过程学到了不少书本上难于获得的知识。它加强了我们的动手能力，分析问题解决问题的能力，同时对我们的学习方法也是一个巨大的锻炼。书本上的知识，需要掌握是比较容易的，但我们所需要的是应用这些知识去解决实际的问题，此次实验给我们提供了一个比较好的锻炼机会。