21摄影测量实习报告

单张相片空间后方交会实习报告一、实习时间实习周数：1周学分数：1分实习时间：第18周二、实习目的1．深入理解单张影像空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。2．利用VisualC++编写一个完整的单张影像空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并进行评定精度。3．通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习内容，培养上机调试程序的动手能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。4、掌握空间后方交会的定义和实现算法(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。(2)算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。5、了解空间后方交会的基本过程(1)获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。(3)确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。(4)计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。(5)逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。(6)逐点计算误差方程式的系数和常数项，组成误差方程式。(7)计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。(8)解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。(9)检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.1′，当3个改正数均小于0.1′时，迭代结束。否则用新的近似值重复(4)～（8）步骤的计算，直到满足要求为止。三、实习环境1．硬件环境：Window操作系统；2．软件环境：VC++四、实习内容利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素并进行精度评定。五、实习原理共线方程式为：)()()()()()()()()()()()(333222333111ssssssssssssZZcYYbXXaZZcYYbXXafyZZcYYbXXaZZcYYbXXafx（1）为了便于计算，将上式线性化，得：()()SSSSSSSSSSSSxxxxxxxxdXdYdZdddXYZyyyyyyyydXdYdZdddXYZ（2）式中，(),()xy为函数的近似值，,,,,,SSSdXdYdZddd为6个外方位元素的改正数。在竖直摄影情况下，角元素都是小角(3)，可用0ASZZH及，将共线方程改写为下面的形式：2222()(1)()(1)SSSSfxxxyxxdXdZfddydHHfffyxyxyydYdZdfdxdHHff（3）当把控制点坐标作为真值，像点坐标作为观测值时，由式（3）列出的误差方程式为：2222(1)(1)xSSxySSyfxxxyVdXdZfddydlHHfffyxyyVdYdZdfdxdlHHff（4）用矩阵形式表示为：VAXl（5）根据最小二乘间接平差原理，按等精度量测，可列出法方程为：TTAAXAL法方程的解为：1()TTXAAAL（6）2.精度评定权倒数为：1()XXTQAA未知数的中误差为：0iiimmQ单位权中误差为：026VVmn式中，代表观测值的点数；6为未知数的个数；2n－6为多余观测数；V为观测点的误差，即像点的观测值与计算值之差。六、程序过程框图：开始获取已知数据：比例尺1/M地面点坐标：像点坐标：x(4),y(4)内方位元素初始值：Xs0,Ys0,Zs0,Phi0,Ome0,Kap0定义旋转矩阵R计算x(i),y(i)的近视值Jinsi(i),Jinsi(i)组成矩阵A(8,6)和L(8,1)计算内方位元素的改正值：dXs,dYs,dZs,dp,do,dk判断内方位元素改正值是否收敛输出Xs,Ys,Zs,φ,ω,κ结束是否实现程序#includeiostream.h#includestdio.h#includestdlib.h#includemath.h#defineN4voidmult(double*m1,double*m2,double*result,inti_1,intj_12,intj_2)//矩阵相乘{inti,j,k;for(i=0;ii_1;i++)for(j=0;jj_2;j++){result[i*j_2+j]=0.0;for(k=0;kj_12;k++)result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];}return;}intinvers_matrix(double*m1,intn)//矩阵求逆{int*is,*js;inti,j,k,l,u,v;doubletemp,max_v;is=(int*)malloc(n*sizeof(int));js=(int*)malloc(n*sizeof(int));if(is==NULL||js==NULL){printf(outofmemory!\n);return(0);}for(k=0;kn;k++){max_v=0.0;for(i=k;in;i++)for(j=k;jn;j++){temp=fabs(m1[i*n+j]);if(tempmax_v){max_v=temp;is[k]=i;js[k]=j;}}if(max_v==0.0){free(is);free(js);printf(inversisnotavailble!\n);return(0);}if(is[k]!=k)for(j=0;jn;j++){u=k*n+j;v=is[k]*n+j;temp=m1[u];m1[u]=m1[v];m1[v]=temp;}if(js[k]!=k)for(i=0;in;i++){u=i*n+k;v=i*n+js[k];temp=m1[u];m1[u]=m1[v];m1[v]=temp;}l=k*n+k;m1[l]=1.0/m1[l];for(j=0;jn;j++)if(j!=k){u=k*n+j;m1[u]*=m1[l];}for(i=0;in;i++)if(i!=k)for(j=0;jn;j++)if(j!=k){u=i*n+j;m1[u]-=m1[i*n+k]*m1[k*n+j];}for(i=0;in;i++)if(i!=k){u=i*n+k;m1[u]*=-m1[l];}}for(k=n-1;k=0;k--){if(js[k]!=k)for(j=0;jn;j++){u=k*n+j;v=js[k]*n+j;temp=m1[u];m1[u]=m1[v];m1[v]=temp;}if(is[k]!=k)for(i=0;in;i++){u=i*n+k;v=i*n+is[k];temp=m1[u];m1[u]=m1[v];m1[v]=temp;}}free(is);free(js);return(1);}voidtranspose(double*m1,double*m2,intm,intn)//矩阵转置{inti,j;for(i=0;im;i++)for(j=0;jn;j++)m2[j*m+i]=m1[i*n+j];return;}voidmain(){doubleXs,Ys,Zs,q,w,k;doublea[3],b[3],c[3];doublex0,y0,f;doublex[N],y[N];doubleX[N],Y[N],Z[N];doublex1[N],y1[N];doublem;doubleL[2*N];doubleXX[6];doubleA[2*N][6];doubleX0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];inti,n=0;doublesum=0,m0;/*---------------输入点地面坐标---------------------*/for(i=0;iN;i++){printf(请输入第%d个点的地面坐标：,i+1);scanf(%lf%lf%lf,&X[i],&Y[i],&Z[i]);}/*---------------输入点像片坐标---------------------*/for(i=0;iN;i++){printf(请输入第%d个点的像片坐标：,i+1);scanf(%lf%lf,&x[i],&y[i]);}coutendl;/*-----------------设定外方位元素初始值--------------*/x0=0;y0=0;f=150;m=30000;Xs=0;Ys=0;Zs=f*m/1000;q=0;w=0;k=0;XX[3]=1;/*------------------迭代计算--------------------------*/while((XX[3]0.00001||XX[4]0.00001||XX[5]0.00001)&&n100){/*----------------旋转矩阵R-----------------------*/a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);/*-----------------像点坐标计算值------------------*/for(i=0;iN;i++){X0[i]=a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs);Y0[i]=a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs);Z0[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);x1[i]=x0-f*X0[i]/Z0[i];y1[i]=y0-f*Y0[i]/Z0[i];}/*-------------误差方程中各偏导数的值--------------*/for(i=0;iN;i++){A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)