1、空间后方交会实习报告遥感信息工程学院 沈翔第一部分:空间后方交会的基本原理空间后方交会的基本原理： 已知影像范围内一定数量的控制点的地面坐标 (或地面摄影坐标)及相对应的影像坐标 ,求的该影像的外方位元素的方法。主要分为：共线条件方程法、角锥体法、直接线性变换法。本程序采用共线条件方程法。第二部分:空间后方交会的计算机程序框图是 否 否否输入原始数据（包括像点坐标观测值、控制点坐标、内方位元素）计算像点坐标 x，y（包括内定向、系统误差改正）计算和确定初值 Xs0, Ys0, Zs0 0, 0, 0、 0, 0, 0组成旋转矩阵 R计算（x） （y）和 lx，ly逐点组成误差方程式并法化迭代次
2、数小于限差否？所有点完否？解法方程，求未知数改正数计算改正后的外方位元素未知数改正数 #include /函数模型void GetR(double phi,double omega,double kappa,double R33);/求 R 阵void Earth2Photo(double earthN3,double R33,double x,double y,double x0,double y0,double f,double xs,double ys,double zs);/计算控制点像点坐标近似值void GetA(double a26,double R33,double earth
3、N3,double photoN2,double x0,double y0,double f,double phi,double omega,double kappa,double zs,int i);/求 A 阵void GetL(double l2,double photoN2,double xN,double yN,int i);/求 L 阵void MatrixT(double a26,double at62,int m,int n);/矩阵转秩double *MatrixMutiple(double *a,int arow,int acol,double *b,int brow,in
4、t bcol);/矩阵相乘double MatrixDet(double x,int n); /求矩阵行列式 double *MatrixInver(double *x,int n);/求逆矩阵/void GetR(double phi,double omega,double kappa,double R33)/求 R 阵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);R1
5、0=cos(omega)*sin(kappa);R11=cos(omega)*cos(kappa);R12=-sin(omega);R20=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);R21=-sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa);R22=cos(phi)*cos(omega);void Earth2Photo(double earthN3,double R33,double x,double y,double x0,double y0,double f,double xs
6、,double ys,double zs)/计算控制点像点坐标近似值int i;for(i=0;ip)z(i-1)*(n-1)+j-1=xi*n+j;if(p%2=0)ss=1;else ss=-1;det+=xp*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 doublen*n;det=MatrixDet(x,n);for(i=0;ij) x_1p*(n-1)+q-
7、1=xp*n+q;if(pi)det_1=MatrixDet(x_1,n-1);if(i+j)%2=0) ss=1;else ss=-1;xxj*n+i=ss*det_1/det;return xx; int main()/定义变量int i,j=0,n=0;/i,j:循环次数/n:收敛次数double photoN2=-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,0.01046,0.06443;/像点坐标double earthN3=36589.41,25273.32,2195.17,37631.08,31324.51,728.69
8、,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 R33;/R 阵元素double xN,yN;/计算的像点坐标double a26,at62;/A 阵及 A 阵的转帙阵double l2;/L 阵double *ata=new double36,*ata_=new double36,*atl
9、=new double6;/指向中间结果的指针 double *X;/改正数的指针double v;/v 的值double m0,mi6;/单位权中误差及未知数的精度/输出初始值cout=delta|*(X+4)=delta|*(X+5)=delta)cout“外方位线元素:tXs(m)ttYs(m)ttZs(m) “endl;cout“tt“xs“ t“ys“ t“zsendl;cout“外方位角元素:tphi(弧度)tomega(弧度)tkappa(弧度)“endl;cout“tt“phi“t“omega“t“kappaendl;GetR(phi,omega,kappa,R);couten
10、dl“R 阵为:“;for(i=0;i3;i+)cout“tt“;for(j=0;j3;j+)coutRij“t“;coutendl;return 0;m0=sqrt(v/(2*N-6);for(j=0;j6;j+)mij=m0*sqrt(*(ata_+7*j);/结果输出cout“= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 求解值= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =“endl;cout“循环次
11、数为:“t“nendl;cout“单位权中误差为:“t“m0endlendl;cout“外方位线元素:tXs(m)ttYs(m)ttZs(m) “endl;cout“tt“xs“ t“ys“ t“zsendl;cout“其精度为:t“;coutmi0“ t“mi1“ t“mi2endlendl;cout“外方位角元素:tphi(弧度)tomega(弧度)tkappa(弧度)“endl;cout“tt“phi“t“omega“t“kappaendl;cout“其精度为:t“;coutmi3“t“mi4“t“mi5endl;GetR(phi,omega,kappa,R);coutendl“R 阵为:“;for(i=0;i3;i+)cout“tt“;for(j=0;j3;j+)coutRij“t“;coutendl;return 1; PRS021 沈翔 200232590166 2004.11.1
