作业4空间后方交会.docx
- 文档编号:6714399
- 上传时间:2023-01-09
- 格式:DOCX
- 页数:21
- 大小:77.63KB
作业4空间后方交会.docx
《作业4空间后方交会.docx》由会员分享,可在线阅读,更多相关《作业4空间后方交会.docx(21页珍藏版)》请在冰豆网上搜索。
作业4空间后方交会
作业报告
空间后方交会
专业:
测绘工程
班级:
2008级
(1)班
**********
*******
2010年4月16日
1作业任务------------------------------------------------------------------------------------3
2作业思想---------------------------------------------------------------------------------------3
3作业条件及数据--------------------------------------------------------------------3
4作业过程---------------------------------------------------------------------------3
5源程序-----------------------------------------------------------------------------4
6计算结果---------------------------------------------------------------------------17
7心得体会与建议-----------------------------------------------------------------------------17
1作业任务
计算近似垂直摄影情况下后方交会解。
即利用摄影测量空间后方交会的方法,获取相片的6个外方位元素。
限差为0.1。
2作业思想
利用摄影测量空间后方交会的方法求解。
该方法的基本思想是利用至少三个一直地面控制点的坐标A(XA,YA,ZA)、B(XB,YB,ZB)C(XC,YC,ZC),与其影像上对应的三个像点的影像坐标a(xa,ya)、b(xb,yb)、c(xc,yc),根据共线方程,反求该相片的外方位元素XS、YS、ZS、φ、ω、κ。
3作业条件及数据
已知摄影机主距f=153.24mm,四对点的像点坐标与相应的地面坐标列入下表:
表1
点号
像点坐标
地面坐标
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
4作业过程
4.1获取已知数据
相片比例尺1/m=1:
10000,内方位元素f=153.24mm,x0,y0;获取控制点的地面测量坐标Xt、Yt、Zt。
4.2量测控制点的像点坐标:
本次作业中为已知。
见表1。
4.3确定未知数的初始值:
在近似垂直摄影情况下,胶原素的初始值为0,即φ0=ω0=κ0=0;线元素中,ZS0=H=mf=1532.4m,XS0、YS0的取值可用四个控制点坐标的平均值,即:
XS0=
=38437.00
YS0=
=89106.62
4.4计算旋转矩阵R:
利用胶原素的近似值计算方向余弦值,组成R阵。
4.5逐点计算像点坐标的近似值:
利用未知数的近似值按共线方程式计算控制点像点坐标的近似值(x)(y)。
4.6组成误差方程:
逐点计算误差方程式的系数和常数项。
4.7组成法方程式:
计算法方程的系数矩阵ATA与常数项ATL。
4.8求解外方位元素:
根据法方程,由式X=(AtA)-1ATL解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。
4.9求解外方位元素:
将求得的外方位元素的改正数与规定的限差(0.1)比较,小于限差则计算终止,否则用新的近似值重复第4.4至4.8步骤的计算,知道满足要求为止。
5源程序
#include
#include
#include
constdoublePRECISION=1e-5;
typedefdoubleDOUBLE[5];
intInputData(int&Num,DOUBLE*&Data,double&m,double&f);
intResection(constint&Num,constDOUBLE*&Data,constdouble&m,constdouble&f);
intInverseMatrix(double*matrix,constint&row);
intmain(intargc,char*argv[])
{
DOUBLE*Data=NULL;
intNum;
doublef(0),m(0);
if(InputData(Num,Data,m,f))
{
if(Data!
=NULL)
{
delete[]Data;
}
return1;
}
if(Resection(Num,Data,m,f))
{
if(Data!
=NULL)
{
delete[]Data;
}
return1;
}
if(Data!
=NULL)
{
delete[]Data;
}
printf("解算完毕...\n");
do{
printf("计算结果保存于\"结果.txt\"文件中\n"
"请选择操作(输入P打开结果数据,R打开原始数据,其它退出程序):
");
fflush(stdin);//刷新输入流
charorder=getchar();
if('P'==order||'p'==order)
{
system("结果.txt");
}
elseif('R'==order||'r'==order)
{
system("data.txt");
}
else
break;
system("cls");
}while
(1);
system("PAUSE");
return0;
}
/**********************************************
*函数名:
InputData
*函数介绍:
从文件(data.txt)中读取数据,
*文件格式如下:
*点数m(未知写作0)
*内方位元素(fx0y0)
*编号xyXYZ
*实例:
40
153.2400
1-86.15-68.9936589.4125273.322195.17
2-53.4082.2137631.0831324.51728.69
3-14.78-76.6339100.9724934.982386.50
410.4664.4340426.5430319.81757.31
*参数:
(in/out)Num(点数),
*(in/out)Data(存放数据),m,f,x0,y0
*返回值:
int,0成功,1文件打开失败,2控制点个
*数不足,3文件格式错误
**********************************************/
intInputData(int&Num,DOUBLE*&Data,double&m,double&f)
{
doublex0,y0;
FILE*fp_input;
if(!
(fp_input=fopen("data.txt","r")))
{
return1;
}
fscanf(fp_input,"%d%lf",&Num,&m);
if(Num<4)
{
return2;
}
fscanf(fp_input,"%lf%lf%lf",&f,&x0,&y0);
f/=1000;
if(m<0||f<0)
{
return3;
}
Data=newDOUBLE[Num];
double*temp=newdouble[Num-1];
doublescale=0;
inti;
for(i=0;i { //读取数据,忽略编号 if(fscanf(fp_input,"%*d%lf%lf%lf%lf%lf", &Data[i][0],&Data[i][1],&Data[i][2], &Data[i][3],&Data[i][4])! =5) { return3; } //单位换算成m Data[i][0]/=1000.0; Data[i][1]/=1000.0; } //如果m未知则归算其值 if(0==m) { for(i=0;i { temp[i]=(Data[i][2]-Data[i+1][2])/(Data[i][0]-Data[i+1][0])+ (Data[i][3]-Data[i+1][3])/(Data[i][1]-Data[i+1][1]); scale+=temp[i]/2.0; } m=scale/(Num-1); } fclose(fp_input); delete[]temp; return0; } /********************************************** *函数名: MatrixMul *函数介绍: 求两个矩阵的积, *参数: Jz1(第一个矩阵),row(第一个矩阵行数), *Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个 *矩阵列数),(out)JgJz(存放结果矩阵) *返回值: void **********************************************/ voidMatrixMul(double*Jz1,constint&row,double*Jz2, constint&line,constint&com,double*JgJz) { for(inti=0;i { for(intj=0;j { doubletemp=0; for(intk=0;k { temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j)); } *(JgJz+i*line+j)=temp; } } } /********************************************** *函数名: OutPut *函数介绍: 向结果.txt文件输出数据 *参数: Q协因数阵,m精度,m0单位权中误差,6个外 *方位元素,旋转矩阵 *返回值: int,0成功,1失败 **********************************************/ intOutPut(constdouble*&Q,constdouble*&m,constdouble&m0, constdouble&Xs,constdouble&Ys,constdouble&Zs, constdouble&Phi,constdouble&Omega, constdouble&Kappa,constdouble*R) { FILE*fp_out; if(! (fp_out=fopen("结果.txt","w"))) { return1; } FILE*fp_input; if(! (fp_input=fopen("data.txt","r"))) { return1; } fprintf(fp_out,"**************************************" "**************************************" "**************************************" "*********************************\n"); fprintf(fp_out,"\n空间后方交会程序(C\\C++)\n测绘一班\n" "学号: 20080729\n姓名: 陈闻亚\n\n"); fprintf(fp_out,"**************************************" "**************************************" "**************************************" "*********************************\n"); fprintf(fp_out,"已知数据: \n\n已知点数: "); intnum;doubletemp,x,y; fscanf(fp_input,"%d%lf",&num,&temp); fprintf(fp_out,"%d\n",num); fprintf(fp_out,"摄影比例尺(0表示其值位置): "); fprintf(fp_out,"%10.0lf\n",temp); fprintf(fp_out,"内方位元素(fx0y0): "); fscanf(fp_input,"%lf%lf%lf",&temp,&x,&y); fprintf(fp_out,"%10lf\t%10lf\t%10lf\n",temp,x,y); for(inti=0;i { doubletemp[5]; fscanf(fp_input,"%*d%lf%lf%lf%lf%lf", &temp[0],&temp[1],&temp[2],&temp[3],&temp[4]); fprintf(fp_out,"%3d\t%10lf\t%10lf\t%10lf\t%10lf\t%10lf\n", i+1,temp[0],temp[1],temp[2],temp[3],temp[4]); } fclose(fp_input); fprintf(fp_out,"**************************************" "**************************************" "**************************************" "*********************************\n"); fprintf(fp_out,"计算结果如下: \n\n外方位元素: \n"); fprintf(fp_out,"\tXs=%10lf\n",Xs); fprintf(fp_out,"\tYs=%10lf\n",Ys); fprintf(fp_out,"\tZs=%10lf\n",Zs); fprintf(fp_out,"\tPhi=%10lf\n",Phi); fprintf(fp_out,"\tOmega=%10lf\n",Omega); fprintf(fp_out,"\tKappa=%10lf\n\n",Kappa); fprintf(fp_out,"旋转矩阵: \n"); for(i=0;i<3;i++) { fprintf(fp_out,"\t"); for(intj=0;j<3;j++) { fprintf(fp_out,"%10lf\t",*(R+i*3+j)); } fprintf(fp_out,"\n"); } fprintf(fp_out,"\n单位权中误差: %10lf\n\n",m0); fprintf(fp_out,"协因数阵: \n"); for(i=0;i<6;i++) { fprintf(fp_out,"\t"); for(intj=0;j<6;j++) { fprintf(fp_out,"%20lf\t",*(Q+i*6+j)); } fprintf(fp_out,"\n"); } fprintf(fp_out,"\n外方位元素精度: "); for(i=0;i<6;i++) { fprintf(fp_out,"%10lf\t",m[i]); } fprintf(fp_out,"\n"); fprintf(fp_out,"**************************************" "**************************************" "**************************************" "*********************************\n"); fclose(fp_out); return0; } /********************************************** *函数名: Resection *函数介绍: 计算 *参数: Num(点数),Data(数据),m,,f(焦距),x0,y0 *返回值: int,0成功,其它失败 **********************************************/ intResection(constint&Num,constDOUBLE*&Data,constdouble&m, constdouble&f) { doubleXs=0,Ys=0,Zs=0; inti,j; //设置初始值 for(i=0;i { Xs+=Data[i][2]; Ys+=Data[i][3]; } Xs/=Num; Ys/=Num; Zs=m*f; doublePhi(0),Omega(0),Kappa(0); doubleR[3][3]={0.0}; double*L=newdouble[2*Num]; typedefdoubleDouble6[6]; Double6*A=newDouble6[2*Num]; double*AT=newdouble[2*Num*6]; double*ATA=newdouble[6*6]; double*ATL=newdouble[6]; double*Xg=newdouble[6]; //迭代计算 do { //旋转矩阵 R[0][0]=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa); R[0][1]=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa); R[0][2]=-sin(Phi)*cos(Omega); R[1][0]=cos(Omega)*sin(Kappa); R[1][1]=cos(Omega)*cos(Kappa); R[1][2]=-sin(Omega); R[2][0]=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa); R[2][1]=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa); R[2][2]=cos(Phi)*cos(Omega); for(i=0;i { doubleX=R[0][0]*(Data[i][2]-Xs)+R[1][0]*(Data[i][3]-Ys)+ R[2][0]*(Data[i][4]-Zs); doubleY=R[0][1]*(Data[i][2]-Xs)+R[1][1]*(Data[i][3]-Ys)+ R[2][1]*(Data[i][4]-Zs); doubleZ=R[0][2]*(Data[i][2]-Xs)+R[1][2]*(Data[i][3]-Ys)+ R[2][2]*(Data[i][4]-Zs); doublexxx,yyy; xxx=-f*X/Z; yyy=-f*Y/Z; //常数项 L[2*i]=Data[i][0]-(-f*X/Z); L[2*i+1]=Data[i][1]-(-f*Y/Z); A[2*i][0]=(R[0][0]*f+R[0][2]*(xxx))/Z; A[2*i][1]=(R[1][0]*f+R[1][2]*(xxx))/Z; A[2*i][2]=(R[2][0]*f+R[2][2]*(xxx))/Z; A[2*i][3]=(yyy)*sin(Omega)-(((xxx)/f)* ((xxx)*cos(Kappa)-(yyy)*sin(Kappa))+ f*cos(Kappa))*cos(Omega); A[2*i][4]=-f*sin(Kappa)-((xxx)/f)*((xxx)* sin(Kappa)+(yyy)*cos(Kappa)); A[2*i][5]=(yyy); A[2*i+1][0]=(R[0][1]*f+R[0][2]*(yyy))/Z; A[2*i+1][1]=(R[1][1]*f+R[1][2]*(yyy))/Z; A[2*i+1][2]=(R[2][1]*f+R[2][2]*(yyy))/Z; A[2*i+1][3]=-(xxx)*sin(Omega)-(((yyy)/f)* ((xxx)*cos(Kappa)-(yyy)*sin(Kappa))- f*sin(Kappa))*cos(Omega); A[2*i+1][4]=-f*cos(Kappa)-((yyy)/f)*((xxx)* sin(Kappa)+(yyy)*cos(Kappa)); A[2*i+1][5]=-(xxx); } //求矩阵A的转置矩阵AT for(i=0;i<2*Num;i++) { for(j=0;j<6;j++) { *(AT+j*2*Num+i)=A[i][j]; } } //求ATA MatrixMul(AT,6,&A[0][0],6,2*Num,ATA); if(InverseMatrix(ATA,6)) return1; MatrixMul(AT,6,L,1,2*Num,ATL); MatrixMul(ATA,6,ATL,1,6,Xg); Xs+=Xg[0]; Ys+=Xg[1]; Zs+=Xg[2]; Phi+=Xg[3]; Omega+=Xg[4]; Kappa+=X
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 作业 空间 后方 交会