单像空间后方交会实验报告c++新版.docx

单像空间后方交会实验报告c++新版.docx

《单像空间后方交会实验报告c++新版.docx》由会员分享，可在线阅读，更多相关《单像空间后方交会实验报告c++新版.docx（20页珍藏版）》请在冰豆网上搜索。

单像空间后方交会实验报告c++新版

单

像

空

间

后

方

交

会

姓名：

学号：

时间：

Echodidthisforyou.2013/4/25

一、作业任务

已知条件：

摄影机主距f=153.24mm，x0=0，y0=0,像片比例尺为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.用三个角元素的初始值按下式计算各方向余弦值，组成旋转矩阵

矩阵中各元素的计算公式如下：

5.逐点计算像点坐标的近似值。

利用未知数的近似值和控制点的地面坐标，带入以下共线方程式，

逐点计算像点坐标的近似值

、

6.逐点计算误差方程式的系数和常数项，组成误差方程式。

由常数项计算公式：

得到常数项矩阵计算式为：

7.计算法方程的是系数矩阵

和常数项

，组成法方程式。

8.解法方程，求得外方位元素的改正数

。

9.用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。

式中：

代表迭代次数。

10.将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。

否则用新的近似值重复4~9满足要求为止。

11.由误差方程的计算式：

以及单位权中误差的计算式

（式中：

表示多余观测数）

以及平差中6个参数的协因数阵

最终得到参数

的中误差为

三、算法流程

四、源程序

源程序代码请见附页

五、计算结果

在经过三次迭代之后得到的最终成果如表1所示

表1，最终计算结果

Xs（M）

Ys（M）

Zs（M）

ψ（弧度）

ω（弧度）

κ（弧度）

39795.452258

27476.462385

7572.685988

-0.003987

0.002114

0.067578

六、结果分析

由计算结果可知在拍摄照片瞬间，摄影中心在地面摄影测量坐标系中的坐标为（39795.452258，27476.462385，7572.685988）（单位：

M），航向倾角ψ为-0.003987弧度，旁向倾角ω为0.002114弧度，像片旋角κ为0.067578弧度。

表2，精度评定结果

参数

Xs

Ys

Zs

ψ

ω

κ

中误差

1.1073

1.2494

0.4881

0.0002

0.0002

0.0001

七、心得与体会

通过本次实验，我对单张像片的空间后方交会的计算原理及实现过程有了很深的认识，并在牢记各种计算公式的推导过程基础上做到了熟练应用。

在本次实验中，我遇到了很多困难，但是都一一得到了解决，在解决困难的过程中的编程能力得到了提高，对其所涉及到的知识的印象也得到了加深。

这个程序本身还有一些不足，例如:

在迭代过程中的矩阵的相乘、输出等，没有采用编写函数的形式，而是逐一进行计算，这是程序整体看起来较臃肿。

在编程过程中为了确保计算的准确性，我还利用matlab同步编写了一个相同功能的计算单张相片空间后方交会的迭代程序，并以此来验算原c++程序的正确性。

现在编写完这个程序，心中很自豪，也有一些遗憾，因为我想用c#来编写，因为c#在暑期实习用过之后，时隔近一年在没有复习，不免把学会的知识又还给了老师，本想通过这次实验对c#进行复习，但是由于时间原因以及一些个人因素，没能达到自己原本的目标。

总之，感谢老师给了我这次锻炼自己、提升自己的机会。

八、附页

1.c++程序

#include"iostream"

usingnamespacestd。

#include"fstream"

#include

constintn=6。

voidinverse（doublec[n][n]）。

templatevoidtranspose（T1*mat1,T2*mat2,inta,intb）。

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

intmain（）

{

doublex[4][2]={-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,0.01046,0.06443}。

doubleX[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:

40000

doubleR[3][3]={0}。

//初始化旋转矩阵R

doubleapprox_x[8]={0}。

//用于存放像点估计值。

doubleA[8][6]={0}。

//定义了一个系数阵

doubleAT[6][8]={0}。

//定义了A的转置矩阵

doubleL[8]={0}。

//定义常数项

constdoublepi=3.14932。

doubleAsum[6][6]={0}。

doubleresult2[6]={0}。

doubleresult1[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"<

}

else

{

cout<<"y"<

}

}

cout<<"已知地面四个点的坐标为：

\n"。

for（i=0。

i<4。

i++）

for（j=0。

j<3。

j++）

{

if（j==0）

{

cout<<"X"。

cout<

cout<<"="<

}

else

if（j==1）

{

cout<<"Y"。

cout<

cout<<"="<

}

else

{

cout<<"Z"。

cout<

cout<<"="。

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][1]=-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[1][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][1]=-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]。

//实现result1[6][8]与l[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<<"进行第"<

\n"。

for（i=0。

i<6。

i++）

{

cout<

f7<

}

cout<

f7.close（）。

getchar（）。

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

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

\n"。

cout<

ofstreamf7（"d:

\\A.txt"）。

f7<

:

fixed。

cout.precision（4）。

for（i=0。

i<6。

i++）

{

cout<

f7<

}

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<

f6<

}

cout<

f6<

}

for（i=0。

i<6。

i++）

cout<

f6.close（）。

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

//

getchar（）。

return0。

}

voidinverse（doublec[n][n]）

{inti,j,h,k。

doublep。

doubleq[n][12]。

for（i=0。

i

i++）//构造高斯矩阵

for（j=0。

j

j++）

q[i][j]=c[i][j]。

for（i=0。

i

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

k++,h++）//消去对角线以下的数据

for（i=k+1。

i

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（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

i++）//将对角线上数据化为1

{

p=1.0/q[i][i]。

for（j=0。

j<12。

j++）

q[i][j]*=p。

}

for（i=0。

i

i++）//提取逆矩阵

for（j=0。

j

j++）

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

}

templatevoidtranspose（T1*mat1,T2*mat2,inta,intb）

{inti,j。

for（i=0。

i

i++）

for（j=0。

j

j++）

mat2[j][i]=mat1[i][j]。

return。

}

templatevoidmulti（T1