空间后方交会

空间后方交会 空间后方交会实习 遥感信息工程学院 沈翔 第一部分:空间后方交会的基本原理 空间后方交会的基本原理: 已知影像范围内一定数量的控制点的地面坐标(或地面摄影坐标)及相对应的影像坐标,求的该影像的外方位元素的方法。主要分为:共线条件方程法、角锥体法、直接线性变换法。本程序采用共线条件方程法。 第二部分:空间后方交会的计算机程序框图 输入原始数据(包括像点坐标观测值、控制点坐标、内方位元素) 计算像点坐标x，y(包括内定向、系统误差改正) X, Y, Z ,, ,, , 计算和确定初值s0s0s0000 、,, ,, , 000 组成旋转矩阵R 计算(x)(y)和lx，ly 是 逐点组成误差方程式并法化 否 迭 代所有点完否, 次 数解法方程，求未知数改正数 小 于 限计算改正后的外方位元素 差 否 否 , 未知数改正数<限差, 否 输出中间结果 求R阵，输出全部计算结果 和出错信息 正常结束 非正常结束 其主要公式如下: 共线方程: aXXbYYcZZ111()()(),，,，,sss xfx,,，0 aXXbYYcZZ333()()(),，,，,sss aXXbYYcZZ222()()(),，,，,sss yfy,,，0 aXXbYYcZZ333()()(),，,，,sss 经线性化后，得到的误差方程式为: ,,,,,,xxxxxxxsss(),,,vxxXYZ,,，,，,，,，,，,，, sss,,,XYZ ,,,,,,yyyyyy(),,,xsssvyyXYZ,,，,，,，,，,，,，, ,,,sssXYZ 组成法方程: TTAPAXAPL, 解方程得: TT,1XAAAL,() 第三部分:试验结果 附录:源程序代码 /*本程序为空间后方交会计算程序(在windows 2003及Visual C++6.0平台中编译通过) 原理:由控制点(至少三个)的地面坐标及相应的影像坐标利用共线方程及最小二乘法来求得 影像外方位元素并确定其精度 其已知变量如下: 内方位元素:x0,y0,f 影象比例尺:m 控制点影象坐标:photo[N][2] 控制点地面坐标:earth[N][3] 其主要中间变量如下: 计算的像点坐标: x[N],y[N] A阵元素: a[2][6] L阵元素: l[2] 待求变量如下: 外方位元素:xs,ys,zs,phi,omega,kappa R阵元素:R[3][3] 单位权中误差:m0 精度:mi[6]*/ const int N=4; #include #include //函数模型 void GetR(double phi,double omega,double kappa,double R[3][3]);//求R阵 void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double x0,double y0,double f,double xs,double ys,double zs);//计算控制点像点 坐标近似值 void GetA(double a[2][6],double R[3][3],double earth[N][3],double photo[N][2],double x0,double y0,double f,double phi,double omega,double kappa,double zs,int i);//求A阵 void GetL(double l[2],double photo[N][2],double x[N],double y[N],int i);//求L阵 void MatrixT(double a[2][6],double at[6][2],int m,int n);//矩阵转秩 double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,int bcol);//矩阵相乘 double MatrixDet(double x[],int n); //求矩阵行列式 double *MatrixInver(double *x,int n);//求逆矩阵 //////////////////////////////////////////////////////////////////////////////////////////////////// void GetR(double phi,double omega,double kappa,double R[3][3])//求R阵 { 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); } void Earth2Photo(double earth[N][3],double R[3][3],double x[],double y[],double x0,double y0,double f,double xs,double ys,double zs)//计算控制点像点 坐标近似值 { int i; for(i=0;ip) z[(i-1)*(n-1)+j-1]=x[i*n+j]; } if(p%2==0) ss=1; else ss=-1; det+=x[p]*ss*MatrixDet(z,n-1); } delete []z; } return det; } double *MatrixInver(double *x,int n)//求逆矩阵 { int i,j,p,q,ss; double det,det_1; double *xx=NULL; xx=new double[n*n]; det=MatrixDet(x,n); for(i=0;ij)) x_1[p*(n-1)+q-1]=x[p*n+q]; if((p>i)&&(qi)&&(q>j)) x_1[(p-1)*(n-1)+q-1]=x[p*n+q]; det_1=MatrixDet(x_1,n-1); if((i+j)%2==0) ss=1; else ss=-1; xx[j*n+i]=ss*det_1/det; } } return xx; } int main() { //定义变量 int i,j=0,n=0;//i,j:循环次数//n:收敛次数 double photo[N][2]={-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663, 0.01046,0.06443};//像点坐标 double earth[N][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}};//地面点坐标 double x0=0,y0=0,f=0.15324;//像片内方位元素 double m=50000;//比例尺 double xs,ys,zs,phi,omega,kappa;//像片外方位元素 double delta;//角元素限值 double R[3][3];//R阵元素 double x[N],y[N];//计算的像点坐标 double a[2][6],at[6][2];//A阵及A阵的转帙阵 double l[2];//L阵 double *ata=new double[36],*ata_=new double[36],*atl=new double[6];//指向中间结果的 指针 double *X;//改正数的指针 double v;//v的值 double m0,mi[6];//单位权中误差及未知数的精度 //输出初始值 cout<<"= = = = = = = = = = = = = = = = = = = = = = = = 初始值= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ="<表

关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf