单像空间后方交会实验报告(c++版)
单像空间后方交会
姓名:
学号:
时间:
目录
一、作业任务 ............................................................................................................... - 3 -
二、计算原理 ............................................................................................................... - 3 -
三、算法流程 ............................................................................................................... - 7 -
四、源程序 ................................................................................................................... - 8 -
五、计算结果 ............................................................................................................... - 8 -
六、结果分析 ............................................................................................................... - 8 -
七、心得与体会 ........................................................................................................... - 8 -
八、附页 ....................................................................................................................... - 8 -
1.c++程序 ........................................................................................................... - 8 -
2.C++程序截图.................................................................................................. - 15 -
3.matlb程序..................................................................................................... - 16 -
一、 作业任务 已知条件:
摄影机主距f=153.24mm ,x0=0,y0=0, 像片比例尺为1:40000,有四对点的像点坐标与相应的地面坐标如下表。
以单像空间后方交会方法,求解该像片的外方位元素。
二、 计算原理
1. 获取已知数据。从航摄资料中查取平均航高与摄影机主距;获取控制点的地面测量坐标并转换为地面摄测坐标。
2. 测量控制点的像点坐标并作系统误差改正。
3. 确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下,按如下方法确定初始值,即
01
,,S S
S
X Y X Y
Z
mf Z n
n
n
=
=
=+
∑∑∑ 0000?ωκ===
式中:m 为摄影比例尺分母;n 为控制点个数
4. 用三个角元素的初始值按下式计算各方向余弦值,组成旋转矩阵R
1231
2
31
2
3a a a R b b b c c c ????=??????
矩阵中各元素的计算公式如下:
12
3123
123cos cos sin sin sin cos sin sin sin cos sin cos cos sin cos cos sin sin cos cos sin sin sin sin cos sin cos cos cos a a a b b b c c c ?κ?ωκ?κ?ωκ?ωωκωκ
ω?κ?ωκ?κ?ωκ?ω
=-??=--??=-?
=??
=??=-??=+?
=-+??=? 5. 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标,带入以下共线方程式,
111333222333()()()
()()()()()()()()()
A S A S A S A S A S A S A S A S A S A S A S A S a X X b Y Y c Z Z x f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z -+-+-?
=-?-+-+-??
-+-+-?=-?-+-+-?
逐点计算像点坐标的近似值()x 、()y
1111111
3131312121211
3131311212122
3232322
()()()()()()()
()()()()()()()()()()()()()()()S S S S S S S S S S S S S S S S S S a X X b Y Y c Z Z x f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z a X X b Y Y c Z Z x f a X X b Y Y c Z Z y -+-+-=--+-+--+-+-=--+-+--+-+-=--+-+-=2222223232321313133333333232323333333314()()()()()()()()()()()()()()()()()()()()(()S S S S S S S S S S S S S S S S S S a X X b Y Y c Z Z f a X X b Y Y c Z Z a X X b Y Y c Z Z x f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z a x f -+-+---+-+--+-+-=--+-+--+-+-=--+-+-=-414143434342424244
343434)()()()()()()()()()()()()S S S S S S S S S S S S X X b Y Y c Z Z a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z ?
?????
???
???
?
???
???
-+-+-??-+-+-?
-+-+-?
=-?-+-+-?
6. 逐点计算误差方程式的系数和常数项,组成误差方程式。
1112131415162122
232425263132333435361414243444546286
51425354555636162
636465664717273747576818283
84
85
86a a a a a a a
a a a a a a a a a a a A a a a a a a A A a a a a a a A a a a a a a A a a a a a a a a a a a a ???
????????????????==?????
?????????????
?
? 由常数项计算公式:
111333222333()()()()()()()()()()()()(1,2,3,4)i S i S i S i
xi i S i S i S i yi i S i S i S i i S i S i S a X X b Y Y c Z Z x f l a X X b Y Y c Z Z L l a X X b Y Y c Z Z y f a X X b Y Y c Z Z i -+-+-?
?+?
?-+-+-???
?==???-+-+-?
??+??-+-+-?
?=
得到常数项矩阵计算式为:
1111111313131212121131313112121223232128134()()()()()()()()()()()()()()()()(S S S S S S S S S S S S S S S S a X X b Y Y c Z Z x f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z a X X b Y Y c Z Z x f a X X b Y Y L L L L L ?-+-+-+-+-+--+-+-+-+-+--+-+-+-+-??????==??????322222222323232131313333333323232333333334)()()()()()()()()()()()()()()()()()()()S S S S S S S S S S S S S S S S S S S S c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z a X X b Y Y c Z Z x f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z x +--+-+-+-+-+--+-+-+-+-+--+-+-+-+-+-1414143434342424244343434()()()()()()()()()()()()S S S S S S S S S S S S a X X b Y Y c Z Z f a X X b Y Y c Z Z a X X b Y Y c Z Z y f a X X b Y Y c Z Z ?
?
??
??
?
??????
?
??
???
?
??
?
??
?
??
?
??
?
??
???
?
-+-+-??
+-+-+-???
?-+-+-??
+-+-+-???
?
7. 计算法方程的是系数矩阵T A A 和常数项T A L ,组成法方程式。 8. 解法方程,求得外方位元素的改正数,,,,,S S S dX dY dZ d d d ?ωκ。
9. 用前次迭代取得的近似值,加本次迭代的改正数,计算外方位元素的新值。
111111,,,,K K K K K K K K K
S S S S S S S S S K K K K K K K K K S
S
S
S
S
S
S
S
S
X X dX Y Y dY Z Z dZ d d d ??
?ωωωκκκ
------=+=+=+=+=+=+
式中:K 代表迭代次数。
10. 将求得的外方位元素改正数与规定的限差比较,若小于限差,则迭代结束。否则用新的近似值重复4~9满足要求为止。 11. 由误差方程的计算式:
1112131415160
212223242526x s s s y s s s v a X a Y a Z a a a x x
v a X a Y a Z a a a y y
?ωκ?ωκ?=?+?+?+?+?+?+-??=?+?+?+?+?+?+-?? 以及单位权中误差的计算式
0σ=
(式中:r 表示多余观测数) 以及平差中6个参数的协因数阵
最终得到参数i X 的中误差为
i
x σσ=
1
()T xx Q A A -=
三、算法流程
四、源程序
源程序代码请见附页
五、计算结果
在经过三次迭代之后得到的最终成果如表1所示
表1,最终计算结果
39795.452258 27476.462385 7572.685988 -0.003987 0.002114 0.067578
六、结果分析
由计算结果可知在拍摄照片瞬间,摄影中心在地面摄影测量坐标系中的坐标为(39795.452258,27476.462385,7572.685988)(单位:米),航向倾角ψ为-0.003987弧度,旁向倾角ω为0.002114弧度,像片旋角κ为0.067578弧度。
表2,精度评定结果
中误差 1.1073 1.2494 0.4881 0.0002 0.0002 0.0001
七、心得与体会
通过本次实验,我对单张像片的空间后方交会的计算原理及实现过程有了很深的认识,并在牢记各种计算公式的推导过程基础上做到了熟练应用。在本次实验中,我遇到了很多困难,但是都一一得到了解决,在解决困难的过程中的编程能力得到了提高,对其所涉及到的知识的印象也得到了加深。
这个程序本身还有一些不足,例如:在迭代过程中的矩阵的相乘、输出等,没有采用编写函数的形式,而是逐一进行计算,这是程序整体看起来较臃肿。
在编程过程中为了确保计算的准确性,我还利用matlab同步编写了一个相同功能的计算单张相片空间后方交会的迭代程序,并以此来验算原c++程序的正确性。
现在编写完这个程序,心中很自豪,也有一些遗憾,因为我想用c#来编写,因为c#在暑期实习用过之后,时隔近一年在没有复习,不免把学会的知识又还给了老师,本想通过这次实验对c#进行复习,但是由于时间原因以及一些个人因素,没能达到自己原本的目标。
总之,感谢老师给了我这次锻炼自己、提升自己的机会。
八、附页
1.c++程序
#include"iostream"
using namespace std;
#include"fstream"
#include
const int n=6;
void inverse(double c[n][n]);
template
template
{
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,4 0426.54,30319.81,757.31};
int i,j,num=0;//num为迭代次数
double X0[6]={0};//设定未知数(XS,YS,ZS,三个角度)初始值
double f=0.15324;//摄影机主距f=153.24mm
double a=1/40000.0;//像片比例尺为1:40000
double R[3][3]={0};//初始化旋转矩阵R
double approx_x[8]={0};//用于存放像点估计值。
double A[8][6]={0};//定义了一个系数阵
double AT[6][8]={0};//定义了A的转置矩阵
double L[8]={0};//定义常数项
const double pi=3.1415926535897932;
double Asum[6][6]={0};
double result2[6]={0};
double result1[6][8]={0};
double sumXYZ[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< } else if(j==1) { cout<<"Y";cout< } else { 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][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]*co s(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])+appr ox_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相乘 int k=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]; } ofstream f7("d:\\A.txt"); f7<< std::fixed; cout<<"进行第"< 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]*20626 5.0)>6); cout<<"\n满足限差条件结束循环,最终结果为:\n"; cout< ofstream f7("d:\\A.txt"); f7<< std::fixed; cout.precision(4); for(i=0;i<6;i++) { cout< } f7.close(); //今 double XG[6][1]; for(i=0;i<6;i++) XG[i][0]=result2[i]; double AXG[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< ofstream f6("d:\\what.txt"); // f6<< std::fixed; for(i=0;i<6;i++) { for(int j=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(); return 0; } void inverse(double c[n][n]) { int i,j,h,k; double p; double q[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]; } template { int i,j; for(i=0;i