空间后方交会编程实习报告一实习目的用程序设计语言（Visual C++或者C语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。

本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。

通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。

二实习内容利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。

三实习数据已知航摄仪的内方位元素：fk =153.24mm，x=y=0.0mm，摄影比例尺为1:50000；四实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。

因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。

可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。

单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，t，w，k。

五实习流程（1）获取已知数据。

从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。

（2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。

（3）确定未知数的初始值。

单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：Z 0s=H=m*f+ΣZi/4； X 0s=Σxi/n; Y 0s=ΣYi/n; t=ω=κ=0;式中：m 为摄影比例尺分母；(4)计算旋转矩阵R 。

利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R 阵。

（5）逐点计算像点坐标的近似值。

利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值（x ）、（y ）。

（6）按下式逐点计算误差方程式的系数和常数项，组成误差方程。

(7)计算法方程的系数矩阵A T A 与常数项A TL ，组成法方程； （8）解求外方位元素。

根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。

（9）检查计算是否收敛。

将所求得的外方位元素的改正数与规定的限差比较，通常对t 、ω、κ、Xs 、Ys 、Zs 的改正数Δt ，Δω，Δκ，ΔXs ，ΔYs ，ΔZs 给予限差，当改正数小于限差时，迭代结束。

否则用新的近似值重复（4）——（8）步骤计算，直到满足要求为止。

（10）空间后方交会的精度估计：按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵（A T A ）-1来解求，此时视像点坐标为等精度不相关观测值。

因为A T A ）-1中第i 个主对角线上的元素Qii 就是法方程式中第i 个未知数的权倒数，若单位权中误差为m0，则第i 个未知数的中误差为： mi=m Qii 0当参加空间后方交会的控制点有n 个时，则单位权中误差可按下式计算：62][-±=n vv m六 主要代码与详解void R(double t,double w,double k,double *a,double *b,double *c) {a[0]=cos(t)*cos(k)-sin(t)*sin(w)*sin(k); a[1]=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k); a[2]=-sin(t)*cos(w); b[0]=cos(w)*sin(k); b[1]=cos(w)*cos(k); b[2]=-sin(w);c[0]=sin(t)*cos(k)+cos(t)*sin(w)*sin(k); c[1]=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k); c[2]=cos(t)*cos(w); }//子函数计算旋转矩阵R 。

利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R 阵。

void main(){int i,m,num;double t,w,k,Xs,Ys,Zs,f; //六个外方位元素与焦距fdouble x[N]={-86.15,-53.40,-14.78,10.46},y[N]={-68.99,82.21,-76.63,64.43};//4个像点坐标doubleX[N]={36589.41,37631.08,39100.97,40426.54},Y[N]={25273.32,31324.51,24934.98,303 19.81},Z[N]={2195.17,728.69,2386.50,757.31};//四个控制点的空间坐标m=50000;f=153.24;// 影像比例尺1/m 与焦距f// 以上主要是已知数据的定义doubleV[6]={0},a[3],b[3],c[3],Xo[N],Yo[N],Zo[N],A[6*M],B[6*M],l[M],C[36],D[6],E[8];for(i=0;i<N;i++){x[i]=x[i]/1000.0;y[i]=y[i]/1000.0;//像点坐标单位换算}// 以下为确定未知数的初始值。

Xs=(X[0]+X[1]+X[2]+X[3])/N;Ys=(Y[0]+Y[1]+Y[2]+Y[3])/N;Zs=m*f+(Z[0]+Z[1]+Z[2]+Z[3])/N;// Xs，Ys，Zs 为摄站点的空间坐标初始值 f=f/1000.0;t=w=k=0.0; //角元素的近似值for(num=1;num<10;num++) //num控制循环迭代次数{R(t,w,k,a,b,c);// 计算旋转矩阵R。

利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。

for(i=0;i<N;i++)//下面的循环用来逐点计算误差方程式的系数和常数项，组成误差方程。

{//以下是用共线条件方程计算控制点像点坐标的近似值（x）、（y）。

Xo[i]=-f*(a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y [i]-Ys)+c[2]*(Z[i]-Zs));Yo[i]=-f*(a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y [i]-Ys)+c[2]*(Z[i]-Zs));Zo[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);//Zo便于后面计算A[12*i+0]=(a[0]*f+a[2]*x[i])/Zo[i];A[12*i+1]=(b[0]*f+b[2]*x[i])/Zo[i];A[12*i+2]=(c[0]*f+c[2]*x[i])/Zo[i];A[12*i+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w); A[12*i+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f; A[12*i+5]=y[i];A[12*i+6]=(a[1]*f+a[2]*y[i])/Zo[i]; A[12*i+7]=(b[1]*f+b[2]*y[i])/Zo[i];A[12*i+8]=(c[1]*f+c[2]*y[i])/Zo[i];A[12*i+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w); A[12*i+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f; A[12*i+11]=-x[i]; l[2*i]=x[i]-Xo[i];l[2*i+1]=y[i]-Yo[i];//计算l }zhuanzhi(A,B,8,6); //转置Axiangchen(B,A,C,6,8,6);//求A TAxiangchen(B,l,D,6,8,1); //求常数项A TLqiuni(C,6);//求A TA 的逆xiangchen(C,D,V,6,6,1);//求改正数VXs+=V[0]; Ys+=V[1]; Zs+=V[2]; t+=V[3]; w+=V[4];k+=V[5]; //结果改正if((fabs(V[0])<=0.00001)&&(fabs(V[1])<=0.00001)&&(fabs(V[2])<=0.00001)&&(fa bs(V[3])<=0.00001)&&(fabs(V[4])<=0.00001)&&(fabs(V[5])<=0.00001))break; //限差控制。

检查计算是否收敛。

将所求得的外方位元素的改正数与规定的限差比较。

当改正数小于限差时，迭代结束。

否则用新的近似值重复计算，直到满足要求为止 。

}//以下是计算单位权中误差 double s=0,m0,v[8];xiangchen(A,V,E,8,6,1); //计算AX for(i=0;i<8;i++) {v[i]=E[i]-l[i];//计算AX-L s+=v[i]*v[i];//计算[vv] } m0=sqrt(s/2);// 单位权中误差按下式计算：62][-±=n vv m ，这里n=4.R(t,w,k,a,b,c); //输出结果printf("像主点的空间坐标为：\n");printf("Xs=%.2f\n",Xs);printf("Ys=%.2f\n",Ys);printf("Zs=%.2f\n",Zs);printf("旋转矩阵R为:\n");for(i=0;i<3;i++)printf(" %.5f",a[i]);printf("\n");for(i=0;i<3;i++)printf(" %.5f",b[i]);printf("\n");for(i=0;i<3;i++)printf(" %.5f",c[i]);printf("\n");printf("单位权中误差为:",m0);cout<<m0<<endl;}七实习结果八实习总结通过本次实习我深刻地理解单片空间后方交会的原理，尤其是共线方程。