单向后方交会实验报告.docx

单向后方交会实验报告.docx

《单向后方交会实验报告.docx》由会员分享，可在线阅读，更多相关《单向后方交会实验报告.docx（15页珍藏版）》请在冰点文库上搜索。

单向后方交会实验报告

Southw^JIaotongUnjiverSity

班级：

测绘一班学号:

20133279

日期：

2015426

、计算原理.

二、算法流程.

三、源程序.

四、计算结果.

五、结果分析.

六、心得体会.

13

13

13

、计算原理

已知条件

摄影机主距f=153.24mm,xO=O,yO=O,像片比例尺为1:

40000,有

四对点的像点坐标与相应的地面坐标如下表。

点号

像点坐标

地面坐标

x（mm）

y（mm）

X（m）

Y（m）

Z（m）

1

-86.15

-68.99

36589.41

25273.32

2195.17

2

-53.40

82.21

37631.08

31324.51

728.69

3

-14.78

-76.63

39100.97

24934.98

2386.50

4

10.46

64.43

40426.54

30319.81

757.31

以单像空间后方交会方法，求解该像片的外方位元素。

二、算法流程

（1）获取已知数据。

从航摄资料中差取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影坐标。

（2）量测控制点的像点坐标并作系统误差改正。

（3）确定未知数的初始值。

在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即

（4）用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵R

KK1.KKK1.KK

d,d,

用共线方程进行空间后方交会的程序框如图所示。

输入原始数据

像点坐标计算，系统误差正

确定外方位因素初始值

组成旋转矩阵R

所有像点完否

是C

解法方程，求外方位元素改正数

计算改正后的外方位元素

外方位元素改正数是否小于限差

输出计算成果，计算并结束

、源程序

#inelude”ostream"usingnamespacestd;

#inelude"fstream"#include

constintn=6;

voidinverse（doublec[n][n]）;

templatevoidtranspose（TI*mat1,T2*mat2,inta,intb）;

templatevtypenameTI,typenameT2>voidmulti（T1*mat1,T2*mat2,T2*result,inta,intb,intc）;intmain（）

{

double

x[4][2]={-0.08615,-0.06899,-0.05340,0.08221，-0.01478，-0.07663,0.01046,0.06443};

double

X[4][3]={36589.41,25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98,

2386.50,40426.54,30319.81,757.31};

inti,j,num=0;//num为迭代次数

doubleX0[6]={0};//设定未知数（XS,YS,ZS,三个角度）初始值

doublef=0.15324;//摄影机主距f=153.24mm

doublea=1/40000.0;//像片比例尺为1:

40000doubleR[3][3]={0};//初始化旋转知阵Rdoubleapprox_x[8]={0};//用于存放像点估计值doubleA[8][6]二{0};//定义了一个系数阵doubleAT[6][8]二{0};//定义了A的转置知阵doubleL[8]={0}:

//定义常数项constdoublepi=3.1415926535897932;doubleAsum[6][6]={0};doubleresult2[6]={0};

doubleresultl[6][8]={0};doublesumXYZ[3]={0};cout.precision（5）;

cout<<”已知像点坐标为:

\n”;for（i=0;i<4;i++）

for（j=0;j<2;j++）

{

cout<

if（j==0）

cout<<”x”<

{

”;cout<

cout<<”Z”;cout<

}

}

cout<

for（j=0;j<3;j++）

for（i=0;i<4;i++）

sumXYZ[j]+=X[i][j];

for（i=0;i<2;i++）

X0[i]=sumXYZ[i]/4;//X0,Y0初始化

X0[i]=1/a*f+sumXYZ[2]/4.0;//对Z0进行初始化

do{

R[0][0]=cos（X0[3]）*cos（X0[5]）-sin（X0[3]）*sin（X0[4]）*sin（X0[5]）;

R[0][I]=-cos（X0[3]）*sin（X0[5]）-sin（X0[3]）*sin（X0[4]）*cos（X0[5]）;

R[0][2]=-sin（X0[3]）*cos（X0[4]）;

R[1][0]=cos（X0[4]）*sin（X0[5]）;

R[I1[1]=cos（X0[4]）*cos（X0[5]）;

R[1][2]=-sin（X0[4]）;

R[2][0]=sin（X0[3]）*cos（X0[5]）+cos（X0[3]）*sin（X0[4]）*sin（X0[5]）;

R[2][l]=-sin（X0[3]）*sin（X0[5]）+cos（X0[3]）*sin（X0[4]）*cos（X0[5]）;

R[2][2]=cos（X0[3]）*cos（X0[4]）;

//第一个像点的估计值，第一个点的坐标存放于X[0][0],X[0][1],X[0][2]

approx_x[0]=-f*（R[0][0]*（X[0][0]-X0[0]）+R[1][0]*（X[0][1]-X0[1]）+R[2][0]*（X[0][2]-X0[2]））/（R[0][2]*（X[0][0]-X0[0]）+R[1][2]*（X[0][1]-X0[1]）+R[2][2]*（X[0][2]-X0[2]））;approx_x[1]=-f*（R[0][1]*（X[0][0]-X0[0]）+R[1][1]*（X[0][1]-X0[1]）+R[2][1]*（X[0][2]-X0[2]））/（R[0][2]*（X[0][0]-X0[0]）+R[1][2]*（X[0][1]-X0[1]）+R[2][2]*（X[0][2]-X0[2]））;//第二个像点的估计值,第2个点的坐标存放于X[1][0],X[1][1],X[1][2]approx_x[2]=-f*（R[0][0]*（X[1][0]-X0[0]）+R[1][0]*（X[1][1]-X0[1]）+R[2][0]*（X[1][2]-X0[2]））/（R[0][2]*（X[1][0]-X0[0]）+R[1][2]*（X[1][1]-X0[1]）+R[2][2]*（X[1][2]-X0[2]））;approx_x[3]=-f*（R[0][1]*（X[1][0]-X0[0]）+R[1][1]*（X[1][1]-X0[1]）+R[2][1]*（X[1][2]-X0[2]））/（R[0][2]*（X[1][0]-X0[0]）+R[1][2]*（X[1][1]-X0[1]）+R[2][2]*（X[1][2]-X0[2]））;

//第三个像点的估计值,第3个点的坐标存放于X[2][0],X[2][1],X[2][2]approx_x[4]=-f*（R[0][0]*（X[2][0]-X0[0]）+R[1][0]*（X[2][1]-X0[1]）+R[2][0]*（X[2][2]-X0[2]））/（R[0][2]*（X[2][0]-X0[0]）+R[1][2]*（X[2][1]-X0[1]）+R[2][2]*（X[2][2]-X0[2]））;approx_x[5]=-f*（R[0][1]*（X[2][0]-X0[0]）+R[1][1]*（X[2][1]-X0[1]）+R[2][1]*（X[2][2]-X0[2]））/（R[0][2]*（X[2][0]-X0[0]）+R[1][2]*（X[2][1]-X0[1]）+R[2][2]*（X[2][2]-X0[2]））;//第四个像点的估计值,第4个点的坐标存放于X[3][0],X[3][1],X[3][2]approx_x[6]=-f*（R[0][0]*（X[3][0]-X0[0]）+R[1][0]*（X[3][1]-X0[1]）+R[2][0]*（X[3][2]-X0[2]））/（R[0][2]*（X[3][0]-X0[0]）+R[1][2]*（X[3][1]-X0[1]）+R[2][2]*（X[3][2]-X0[2]））;approx_x[7]=-f*（R[0][1]*（X[3][0]-X0[0]）+R[1][1]*（X[3][1]-X0[1]）+R[2][1]*（X[3][2]-X0[2]））/（R[0][2]*（X[3][0]-X0[0]）+R[1][2]*（X[3][1]-X0[1]）+R[2][2]*（X[3][2]-X0[2]））;for（i=0;i<4;i++）

{

//第i个像点估计值放在approx_x[2*（i-1）],approx_x[2*i-1]/*a11*/A[2*i][0]=（R[0][0]*f+R[0][2]*approx_x[2*i]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;/*a12*/A[2*i][1]=（R[1][0]*f+R[1][2]*approx_x[2*i]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;

/*a13*/A[2*i][2]=（R[2][0]*f+R[2][2]*approx_x[2*i]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;

/*a21*/A[2*i+1][0]=（R[0][1]*f+R[0][2]*approx_x[2*i+1]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;

/*a22*/A[2*i+1][1]=（R[1][1]*f+R[1][2]*approx_x[2*i+1]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;

/*a23*/A[2*i+1][2]=（R[2][1]*f+R[2][2]*approx_x[2*i+1]）/（R[0][2]*（X[i][0]-X0[0]）+R[1][2]*（X[i][1]-X0[1]）+R[2][2]*（X[i][2]-X0[2]））;

/*a14*/A[2*i][3]=approx_x[2*i+1]*sin（X0[4]）-（approx_x[2*i]/f*（approx_x[2*i]*cos（X0[5]）-approx_x[2*i+1]*sin（X0[5]））+f*cos（X0[5]））*cos（X0[4]）;

/*a15*/A[2*i][4]=-f*sin（X0[5]）-approx_x[2*i]/f*（approx_x[2*i]*sin（X0[5]）+approx_x[2*i+1]*cos（X0[5]））;/*a16*/A[2*i][5]=approx_x[2*i+1];

/*a24*/A[2*i+1][3]=-1*approx_x[2*i]*sin（X0[4]）-（approx_x[2*i+1]/f*（approx_x[2*i]*cos（X0[5]）-approx_x[2*i+1]*sin（X0[5]））-f*sin（X0[5]））*cos（X0[4]）;

/*a25*/A[2*i+1][4]=-1*f*cos（X0[5]）-approx_x[2*i+1]/f*（approx_x[2*i]*sin（X0[5]）+approx_x[2*i+1]*cos（X0[5]））;

/*a26*/A[2*i+1][5]=-approx_x[2*i];

}

//进行常数项的初始化

for（i=0;i<4;i++）

{

L[2*i]=x[i][0]-approx_x[2*i];

L[2*i+1]=x[i][1]-approx_x[2*i+1];

}

//A的转置矩阵

for（i=0;i<8;i++）

for（j=0;j<6;j++）

{

AT[j][i]=A[i][j];

}

//实现A与AT相乘

intk=0;

for（i=0;i<6;i++）

for（j=0;j<6;j++）

Asum[i][j]=0;

for（i=0;i<6;i++）

for（k=0;k<6;k++）

for（j=0;j<8;j++）

Asum[i][k]+=AT[i][j]*A[j][k];

//得到AT*A的逆矩阵存放在inverseAsum[6][6]中

inverse（Asum）;

//实现矩阵Asum[6][6]与AT[6][8]的相乘,结果存放在result1[6][8]中for（i=0;i<6;i++）

for（j=0;j<8;j++）

result1[i][j]=0;

for（i=0;i<6;i++）for（k=0;k<8;k++）

for（j=0;j<6;j++）

result1[i][k]+=Asum[i][j]*AT[j][k];

//实现resultl[6][8]与1[8]的相乘，得到结果放在result2[6]中;

for（i=0;i<6;i++）

result2[i]=0;

for（i=0;i<6;i++）

for（j=0;j<8;j++）

result2[i]+=result1[i][j]*L[j];

num++;

for（i=0;i<6;i++）

{

X0[i]=X0[i]+result2[i];

}

ofstreamf7（"d:

\\A.txt"）;

f7<

:

fixed;

cout<<"进行第"<

for（i=0;i<6;i++）

{cout<

}

cout<

f7.close（）;

getchar（）;

}

while（abs（result2[3]*206265.0）>6||abs（result2[4]*206265.0）>6||abs（result2[5]*20626

5.0）>6）;

cout<<"\n满足限差条件结束循环，最终结果为：

\n";

cout<

ofstreamf7（"d:

\\A.txt"）;

f7<

:

fixed;

cout.precision（4）;

for（i=0;i<6;i++）

{

cout<

}

f7.close（）;

//今

doubleXG[6][1];

for（i=0;i<6;i++）

XG[i][0]=result2[i];

doubleAXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];

multi（A,XG,AXG,8,6,1）;

for（i=0;i<8;i++）

//计算改正数

V[i][0]=AXG[i][0]-L[i];

transpose（V,VT,1,8）;

multi（VT,V,VTV,1,8,1）;

m0=VTV[0][0]/2;

cout<

ofstreamf6（"d:

\\what.txt"）;

//f6<

:

fixed;

for（i=0;i<6;i++）

{

for（intj=0;j<6;j++）

{

D[i][j]=m0*Asum[i][j];cout<

}

cout<

f6<

}

for（i=0;i<6;i++）

cout<

//屏幕输出误差方程系数阵、常数项、改正数//

getchar（）;

return0;

}

voidinverse（doublec[n][n]）

{inti,j,h,k;doublep;

doubleq[n][12];

for（i=0;i

for（j=0;j

q[i][j]=c[i][j];

for（i=0;i

for（j=n;j<12;j++）

{

if（i+6==j）

q[i][j]=1;

else

q[i][j]=0;

}

for（h=k=0;k

{if（q[i][h]==0）continue;p=q[k][h]/q[i][h];

for（j=0;j<12;j++）

{

q[i][j]*=p;

q[i][j]-=q[k][j];

}

}

for（h=k=n-1;k>0;k--,h--）//消去对角线以上的数据for（i=k-1;i>=0;i--）

{

if（q[i][h]==0）

continue;

p=q[k][h]/q[i][h];

for（j=0;j<12;j++）

{

q[i][j]*=p;q[i][j]-=q[k][j];

}

}

for（i=0;i

{p=1.0/q[i][i];

for（j=0;j<12;j++）

q[i][j]*=p;

}

for（i=0;i

for（j=0;j

c[i][j]=q[i][j+6];

}

templatevoidtranspose（T1*mat1,T2*mat2,int

a,intb）

{inti,j;

for（i=0;i

for（j=0;j

mat2[j][i]=mat1[i][j];

return;

}

templatevoidmulti（T1*mat1,T2*mat2,T2*result,inta,intb,intc）

{

for（i=0;i

for（k=0;k

return;

程序截图

i.TF農2坎这代爭待到肘占•空.叩，dK灿_教分別.为:

-47.80370-S?

.30563-2Q.S7G4C0■丽046

中i天着为,

I,1»73i.2吐Murna,0&9s

乩加'叶I

四、计算结果

Xs

Ys

Zs

39795.452258

27476.462385

7572.685988

-0.003987

0.002114

0.067578

五、结果分析

参数

Xs

Ys

Zs

中误差

1.1073

1.2494

0.4881

0.0002

0.0002

0.0001

由计算结果可知：

摄影中心在地面摄影测量坐标系中的坐标为航向倾角为-0.003987弧度，旁向倾角

（39795.452258.27476.462385.7572.685988）m

为0.002114弧度，像片旋角为0.067578弧度。

六、心得体会

本次实验是关于单张像片的空间后方交会的计算原理及实现过程的课程作业，其中解题

步骤和公式是基本，熟练运用公式并且成功编程是重点。

很遗憾，编程方面我是弱项，不得

不借鉴网上的编程，至今对编程的某些地方还没弄明白，我将继续学习和熟练编程的方法的

技巧。