作业4空间后方交会

1、摄影测量作业报告之空间后方交会陈闻亚 20080729作业报告空间后方交会专 业： 测绘工程 班 级： 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

2、，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.9936589.4125273.322195.172-53.40-82.2137631.0831324.51.728.69314.7876.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31

3、4作业过程 41 获取已知数据 相片比例尺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 逐点计算像点坐标的近似值： 利用未知数的近似值按共线方程式计算控制点像

4、点坐标的近似值（x）（y）。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 d

5、ouble5;int inputdata(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

6、(data!=null) delete 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=

7、order | r=order) 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

8、.97 24934.98 2386.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,&n

9、um,&m);if (num4) 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 data

10、i0/=1000.0; datai1/=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(第二个矩阵列数

11、),com(第一个*矩阵列数),(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

12、*函数介绍：向结果.txt文件输出数据*参数：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(结果.tx

13、t,w) return 1;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

14、_out,摄影比例尺（0表示其值位置）：);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%10lf

15、t%10lft%10lft%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);fprin

16、tf(fp_out,tkappa=%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

17、+i*6+j); fprintf(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

18、 &m, const double &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

19、;double *atl=new 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(

20、phi)*sin(omega)*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-

21、zs); double xxx,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

22、)* sin(kappa)+(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=

23、-(xxx); /求矩阵a的转置矩阵at for (i=0;i2*num;i+) for (j=0;j=precision |fabs(xg1)=precision | fabs(xg2)=precision |fabs(xg3)=precision | fabs(xg4)=precision | (xg5)=precision);/注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值，/由于变换很小忽略double *q=ata;double *v=new double2*num;matrixmul(&a00,2*num,xg,1,6,v);double vtv=0;for(i=0;i2*

24、num;i+) vi-=li; vtv+=vi*vi;double m0=sqrt(vtv/(2*num-6);double *mm=new double6;for (i=0;i6;i+) mmi=sqrt(*(q+i*6+i)*m0;output(q,mm,m0,xs,ys,zs,phi,omega,kappa,&r00);delete l;delete a;delete at;delete ata;delete atl;delete xg;delete mm;delete v;return 0;void swap(double &a,double &b)double temp=a;a=b;

25、b=temp;/*函数名：inversematrix *函数介绍：求矩阵的逆（高斯-约当法） *输入参数：(in/out)matrix（矩阵首地址），*（in）row（矩阵阶数）*输出参数：matrix（原矩阵的逆矩阵）*返回值：int ，0成功，1失败*调用函数：swap(double&,double&)*/int inversematrix(double *matrix,const int &row)double *m=new doublerow*row;double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;irow;i+) for (j=0;

26、jrow;j+) *pt=*ptemp; ptemp+; pt+; int k;int *is=new introw,*js=new introw;for (k=0;krow;k+) double max=0; /全选主元 /寻找最大元素 for (i=k;irow;i+) for (j=k;jmax) max=*(m+i*row+j); isk=i; jsk=j; if (0 = max) return 1; /行交换 if (isk!=k) for (i=0;irow;i+) swap(*(m+k*row+i),*(m+isk*row+i); /列交换 if (jsk!=k) for (i=0;irow;i+) swap(*(m+i*row+k),*(m+i*row+jsk); *(m+k*row+k)=1/(*(m+k*row+k); for (j=0;jrow;j+) if (j!=k) *(m+k*row+j)*=*(m+k*row+k); for (i=0;irow;i+) if (i!=k) for (j=0;jrow;j+) if(j!=k) *(m+i*row+j)-=*(m+i*row+k)*(m+