1、作业4空间后方交会作业报告空间后方交会专 业: 测绘工程 班 级: 2008级(1)班 姓 名: 陈闻亚 指导教师: 陈强 2010 年 4 月 16 日1 作业任务 - 32 作业思想 - 33 作业条件及数据- 34 作业过程 - 35 源程序 - 46 计算结果 - 177心得体会与建议- 171 作业任务计算近似垂直摄影情况下后方交会解。即利用摄影测量空间后方交会的方法,获取相片的6个外方位元素。限差为0.1。2作业思想利用摄影测量空间后方交会的方法求解。该方法的基本思想是利用至少三个一直地面控制点的坐标A(XA,YA,ZA)、B(XB,YB,ZB)C(XC,YC,ZC),与其影像上对
2、应的三个像点的影像坐标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.9936589.4125273.322195.172-53.40-82.2137631.0831324.51.728.69314.7876.6339100.9724934.982386.50410.4664.4340426.5430319.81757.314作业过程 41 获取已知数
3、据 相片比例尺1/m=1:10000,内方位元素f=153.24mm,x0,y0;获取控制点的地面测量坐标Xt、Yt、Zt。42 量测控制点的像点坐标: 本次作业中为已知。见表1。43 确定未知数的初始值: 在近似垂直摄影情况下,胶原素的初始值为0,即0 = 0 = 0=0;线元素中,ZS0=H=mf=1532.4m,XS0、YS0的取值可用四个控制点坐标的平均值,即:XS0= =38437.00YS0= =89106.6244 计算旋转矩阵R: 利用胶原素的近似值计算方向余弦值,组成R阵。45 逐点计算像点坐标的近似值: 利用未知数的近似值按共线方程式计算控制点像点坐标的近似值(x)(y)。
4、46 组成误差方程: 逐点计算误差方程式的系数和常数项。47 组成法方程式: 计算法方程的系数矩阵ATA与常数项ATL。48 求解外方位元素: 根据法方程,由式X=(AtA)-1 ATL解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。49 求解外方位元素: 将求得的外方位元素的改正数与规定的限差(0.1)比较,小于限差则计算终止,否则用新的近似值重复第4.4至4.8步骤的计算,知道满足要求为止。5 源程序#include #include #include const double PRECISION=1e-5;typedef double DOUBLE5;int Inp
5、utData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f);int InverseMatrix(double *matrix,const int &row);int main(int argc, char* argv)DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f) if (Data!=NULL) d
6、elete Data; return 1;if(Resection(Num,Data,m,f) if (Data!=NULL) delete Data; return 1;if (Data!=NULL) delete Data;printf(解算完毕.n);do printf(计算结果保存于结果.txt文件中n 请选择操作(输入P打开结果数据,R打开原始数据,其它退出程序):); fflush(stdin); /刷新输入流 char order=getchar(); if (P=order | p=order) system(结果.txt); else if (R=order | r=orde
7、r) system(data.txt); else break; system(cls);while(1);system(PAUSE);return 0;/*函数名:InputData *函数介绍:从文件(data.txt)中读取数据,*文件格式如下: *点数 m(未知写作0)* 内方位元素(f x0 y0)*编号 x y X Y Z*实例:4 0153.24 0 01 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2
8、386.504 10.46 64.43 40426.54 30319.81 757.31*参数:(in/out)Num(点数),*(in/out)Data(存放数据),m,f,x0,y0*返回值:int ,0成功,1文件打开失败,2控制点个*数不足,3文件格式错误*/int InputData(int &Num, DOUBLE *&Data,double &m,double &f)double x0,y0;FILE *fp_input;if (!(fp_input=fopen(data.txt,r) return 1;fscanf(fp_input,%d%lf,&Num,&m);if (Num
9、4) return 2;fscanf(fp_input,%lf%lf%lf,&f,&x0,&y0);f/=1000;if (m0 | f0) return 3;Data=new DOUBLENum;double *temp= new doubleNum-1;double scale=0;int i;for (i=0;iNum;i+) /读取数据,忽略编号 if(fscanf(fp_input,%*d%lf%lf%lf%lf%lf, &Datai0,&Datai1,&Datai2, &Datai3,&Datai4)!=5) return 3; /单位换算成m Datai0/=1000.0; Da
10、tai1/=1000.0;/如果m未知则归算其值if (0=m) for (i=0;iNum-1;i+) tempi=(Datai2-Datai+12)/(Datai0-Datai+10)+ (Datai3-Datai+13)/(Datai1-Datai+11); scale+=tempi/2.0; m=scale/(Num-1);fclose(fp_input);delete temp;return 0;/*函数名:MatrixMul *函数介绍:求两个矩阵的积,*参数:Jz1(第一个矩阵),row(第一个矩阵行数),*Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个*矩阵列数
11、),(out)JgJz(存放结果矩阵)*返回值:void*/void MatrixMul(double *Jz1,const int &row,double *Jz2, const int &line,const int &com,double *JgJz)for (int i=0;irow;i+) for (int j=0;jline;j+) double temp=0; for (int k=0;kcom;k+) temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j); *(JgJz+i*line+j)=temp; /*函数名:OutPut *函数介绍:向结果.txt文
12、件输出数据*参数:Q协因数阵,m精度,m0单位权中误差,6个外*方位元素,旋转矩阵*返回值:int,0成功,1失败*/int OutPut(const double *&Q,const double *&m,const double &m0, const double &Xs,const double &Ys,const double &Zs, const double &Phi,const double &Omega, const double &Kappa,const double *R)FILE *fp_out;if (!(fp_out=fopen(结果.txt,w) return 1;
13、FILE *fp_input;if (!(fp_input=fopen(data.txt,r) return 1;fprintf(fp_out,* * * *n);fprintf(fp_out,n空间后方交会程序(CC+)n测绘一班n 学号:20080729n姓名:陈闻亚nn);fprintf(fp_out,* * * *n);fprintf(fp_out,已知数据:nn已知点数:);int num;double temp,x,y;fscanf(fp_input,%d%lf,&num,&temp);fprintf(fp_out,%dn,num);fprintf(fp_out,摄影比例尺(0表示
14、其值位置):);fprintf(fp_out,%10.0lfn,temp);fprintf(fp_out,内方位元素(f x0 y0):);fscanf(fp_input,%lf%lf%lf,&temp,&x,&y);fprintf(fp_out,%10lft%10lft%10lfn,temp,x,y);for (int i=0;inum;i+) double temp5; fscanf(fp_input,%*d%lf%lf%lf%lf%lf, &temp0,&temp1,&temp2,&temp3,&temp4); fprintf(fp_out,%3dt%10lft%10lft%10lft%
15、10lft%10lfn, i+1,temp0,temp1,temp2,temp3,temp4);fclose(fp_input);fprintf(fp_out,* * * *n);fprintf(fp_out,计算结果如下:nn外方位元素:n);fprintf(fp_out,tXs=%10lfn,Xs);fprintf(fp_out,tYs=%10lfn,Ys);fprintf(fp_out,tZs=%10lfn,Zs);fprintf(fp_out,tPhi=%10lfn,Phi);fprintf(fp_out,tOmega=%10lfn,Omega);fprintf(fp_out,tKap
16、pa=%10lfnn,Kappa);fprintf(fp_out,旋转矩阵:n);for (i=0;i3;i+) fprintf(fp_out,t); for (int j=0;j3;j+) fprintf(fp_out,%10lft,*(R+i*3+j); fprintf(fp_out,n);fprintf(fp_out,n单位权中误差:%10lfnn,m0);fprintf(fp_out,协因数阵:n);for (i=0;i6;i+) fprintf(fp_out,t); for (int j=0;j6;j+) fprintf(fp_out,%20lft,*(Q+i*6+j); fprin
17、tf(fp_out,n);fprintf(fp_out,n外方位元素精度:);for (i=0;i6;i+) fprintf(fp_out,%10lft,mi);fprintf(fp_out,n);fprintf(fp_out,* * * *n);fclose(fp_out);return 0;/*函数名:Resection *函数介绍:计算*参数:Num(点数),Data(数据),m,f(焦距),x0,y0*返回值:int,0成功,其它失败*/int Resection(const int &Num,const DOUBLE *&Data,const double &m, const dou
18、ble &f)double Xs=0,Ys=0,Zs=0;int i,j;/设置初始值for (i=0;iNum;i+) Xs+=Datai2; Ys+=Datai3;Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);double R33=0.0;double *L=new double2*Num;typedef double Double66;Double6 *A=new Double62*Num;double *AT=new double2*Num*6;double *ATA=new double6*6;double *ATL=n
19、ew double6;double *Xg=new double6;/迭代计算do /旋转矩阵 R00=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa); R01=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa); R02=-sin(Phi)*cos(Omega); R10=cos(Omega)*sin(Kappa); R11=cos(Omega)*cos(Kappa); R12=-sin(Omega); R20=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega
20、)*sin(Kappa); R21=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa); R22=cos(Phi)*cos(Omega); for (i=0;iNum;i+) double X=R00*(Datai2-Xs)+R10*(Datai3-Ys)+ R20*(Datai4-Zs); double Y=R01*(Datai2-Xs)+R11*(Datai3-Ys)+ R21*(Datai4-Zs); double Z=R02*(Datai2-Xs)+R12*(Datai3-Ys)+ R22*(Datai4-Zs); double xx
21、x,yyy; xxx=-f*X/Z; yyy=-f*Y/Z; /常数项 L2*i=Datai0-(-f*X/Z); L2*i+1=Datai1-(-f*Y/Z); A2*i0=(R00*f+R02*(xxx)/Z; A2*i1=(R10*f+R12*(xxx)/Z; A2*i2=(R20*f+R22*(xxx)/Z; A2*i3=(yyy)*sin(Omega)-(xxx)/f)* (xxx)*cos(Kappa)-(yyy)*sin(Kappa)+ f*cos(Kappa)*cos(Omega); A2*i4=-f*sin(Kappa)-(xxx)/f)*(xxx)* sin(Kappa)+
22、(yyy)*cos(Kappa); A2*i5=(yyy); A2*i+10=(R01*f+R02*(yyy)/Z; A2*i+11=(R11*f+R12*(yyy)/Z; A2*i+12=(R21*f+R22*(yyy)/Z; A2*i+13=-(xxx)*sin(Omega)-(yyy)/f)* (xxx)*cos(Kappa)-(yyy)*sin(Kappa)- f*sin(Kappa)*cos(Omega); A2*i+14=-f*cos(Kappa)-(yyy)/f)*(xxx)* sin(Kappa)+(yyy)*cos(Kappa); A2*i+15=-(xxx); /求矩阵A的转置矩阵AT for (i=0;i2*Num;i+) for (j=0;j6;j+) *(AT+j*2*Num+i)=Aij; /求ATA MatrixMul(AT,6,&A00,6,2*Num,ATA); if(InverseMatrix(ATA,6) return 1; MatrixMul(AT,6,L,1,2*Num,ATL); MatrixMul(ATA,6,ATL,1,6,Xg); Xs+=Xg0; Ys+=Xg1; Zs+=Xg2; Phi+=Xg3; Omega+=Xg4; Kapp
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1