单像空间后方交会实习报告

1摄影测量学单像空间后方交会实习报告2一、实习目的1.掌握空间后方交会的定义和实现算法(1)定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，φ，ω，κ。(2)算法：由于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数△Xs，△Ys，△Zs，△φ，△ω，△κ。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。2.了解摄影测量平差的基本过程(1)获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高）、内方位元素x0，y0，f；获取控制点的空间坐标X，Y，Z。(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。(3)确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。(4)计算旋转矩阵R。利用角元素近似值计算方向余弦值，组成R阵。(5)逐点计算像点坐标的近似值。利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)。(6)逐点计算误差方程式的系数和常数项，组成误差方程式。(7)计算法方程的系数矩阵ATA与常数项ATL，组成法方程式。(8)解求外方位元素。根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。(9)检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.000001弧度，当3个改正数均小于0.000001弧度时，迭代结束。否则用新的近似值重复(4)～（8）步骤的计算，直到满足要求为止。3.通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。4.实习过程：4.1学习单张像片空间后方交会的基本理论，掌握其基本思想。如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。而单像空间后方交会就是用于测定像片的外方位元素的，它的基本思想是：以单幅影像为基础，从影像所覆盖的地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,p,w,k.由于空间后方交会所采用的数学模型共线方程是非线性函数，为了便于外方位元素的解求，首先将其线性化。4.2在纸上绘出空间后方交会的计算机程序框图。为了能够在宏观上指导我们编写程序，我们需要在草稿纸上绘出程序框图。34输入原始数据归算像点坐标x，y计算和确定初值Xs0，Ys0，Zs0，φ0，ω0，κ0组成旋转矩阵R计算(x),(y)和lx，ly逐点组成误差方程式并法化所有点完否？计算改正后的外方位元素解法方程，求未知数改正数未知数改正数限差否？整理并输出计算结果正常结束否迭代次数小于限差否？否是否输出中间结果和出错信息非正常结束5二、源代码usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;usingSystem.Text;namespace单像空间的后方交会{publicclasscalculate{privatedoublej,k,l,Xs,Ys,Zs;//六个外方位元素privatedoublef=28.1539;//主距structpoint//像点和地面点坐标{publicdoublex,y,X,Y,Z;}privatepoint[]p=newpoint[4];//存取控制点的坐标privatedouble[]R=newdouble[9];//旋转矩阵privatedouble[]a=newdouble[8];//近似值坐标privatedouble[]L=newdouble[8];//误差方程常数项privatedouble[,]A=newdouble[8,6];//误差方程系数项privateintcount=0;publicvoidY(double[]q){for(intn=0;n4;n++){intm=n*5;Console.WriteLine(请输入第｛0｝控制点的坐标,n+1);q[m]=Convert.ToDouble(Console.ReadLine());q[m+1]=Convert.ToDouble(Console.ReadLine());q[m+2]=Convert.ToDouble(Console.ReadLine());q[m+3]=Convert.ToDouble(Console.ReadLine());q[m+4]=Convert.ToDouble(Console.ReadLine());p[n].x=q[m];p[n].y=q[m+1];p[n].X=q[m+2];p[n].Y=q[m+3];p[n].Z=q[m+4];}doubleave=0,sum=0;//求比例尺分母，用来求外方位元素for(intn=0;n3;n++){for(intm=n+1;m4;m++)6{sum+=Math.Sqrt(Math.Pow(p[n].X-p[m].X,2)+Math.Pow(p[n].Y-p[m].Y,2))/Math.Sqrt(Math.Pow(p[n].x-p[m].x,2)+Math.Pow(p[n].y-p[m].y,2));}}ave=sum/6;doublej=0.054882;k=0.057034;l=-0.036175;//六个外方位元素的近似值doubleXs=500215.49;doubleYs=4185301.89;doubleZs=1475.56;}privatedoublesin(doublem){returnMath.Sin(m);}privatedoublecos(doublem){returnMath.Cos(m);}publicvoidcalX()//计算旋转矩阵{R[0]=cos(j)*cos(k)-sin(j)*sin(l)*sin(k);R[1]=-cos(j)*sin(k)-sin(j)*sin(l)*cos(k);R[2]=-sin(j)*cos(l);R[3]=cos(l)*sin(k);R[4]=cos(l)*cos(k);R[5]=-sin(l);R[6]=sin(j)*cos(k)+cos(j)*sin(l)*sin(k);R[7]=-sin(j)*sin(k)+cos(j)*sin(l)*cos(k);R[8]=cos(j)*cos(l);}publicvoidcalJSZ()//计算像点坐标近似值{for(intn=0;n4;n++){a[2*n]=-f*(R[0]*(p[n].X-Xs)+R[3]*(p[n].Y-Ys)+R[6]*(p[n].Z-Zs))/(R[2]*(p[n].X-Xs)+R[5]*(p[n].Y-Ys)+R[8]*(p[n].Z-Zs));a[2*n+1]=-f*(R[1]*(p[n].X-Xs)+R[4]*(p[n].Y-Ys)+R[7]*(p[n].Z-Zs))/(R[2]*(p[n].X-Xs)+R[5]*(p[n].Y-Ys)+R[8]*(p[n].Z-Zs));}}publicvoidcalXSJZandCSX()//计算系数矩阵和常数项{7for(intn=0;n4;n++)//计算常数项{L[2*n]=p[n].x-a[2*n];L[2*n+1]=p[n].y-a[2*n+1];}for(intn=0;n4;n++)//计算系数矩阵{doublez=R[2]*(p[n].X-Xs)+R[5]*(p[n].Y-Ys)+R[8]*(p[n].Z-Zs);A[2*n,0]=(R[0]*f+R[2]*p[n].x)/z;A[2*n,1]=(R[3]*f+R[5]*p[n].x)/z;A[2*n,2]=(R[6]*f+R[8]*p[n].x)/z;A[2*n,3]=p[n].y*sin(l)-(p[n].x*(p[n].x*cos(k)-p[n].y*sin(k))/f+f*cos(k))*cos(l);A[2*n,4]=-f*sin(k)-p[n].x/f*(p[n].x*sin(k)+p[n].y*cos(k));A[2*n,5]=p[n].y;A[2*n+1,0]=(R[1]*f+R[2]*p[n].y)/z;A[2*n+1,1]=(R[4]*f+R[5]*p[n].y)/z;A[2*n+1,2]=(R[7]*f+R[8]*p[n].y)/z;A[2*n+1,3]=p[n].x*sin(l)-(p[n].y*(p[n].x*cos(k)-p[n].y*sin(k))/f-f*sin(k))*cos(l);A[2*n+1,4]=-f*cos(k)-p[n].y/f*(p[n].x*sin(k)+p[n].y*cos(k));A[2*n,5]=-p[n].x;}}publicdoublecalJSGZS()//计算改正数{double[,]AT=newdouble[6,8];//A转置double[,]temp=newdouble[6,6];//A转置与A的乘积double[]X=newdouble[6];//改正数double[]ATL=newdouble[6];//A转置与L相乘的积intn,m,s;for(n=0;n8;n++)//求A的转置矩阵ATfor(m=0;m6;m++){AT[m,n]=AT[n,m];}for(n=0;n6;n++)//求A转置与A的乘积for(m=0;m6;m++){temp[n,m]=0;for(s=0;s8;s++){temp[n,m]+=AT[n,s]*A[s,m];}}8MA(temp);//求逆for(n=0;n6;n++)//求A的转置与L的乘积{for(m=0;m8;m++){ATL[n]+=AT[n,m]*L[m];}}for(n=0;n6;n++)//求改正数X{for(m=0;m6;m++){X[n]+=temp[n,m]*ATL[m];}}Xs+=X[0];Ys+=X[1];Zs+=X[2];j+=X[3];k+=X[4];l+=X[5];returnXs;returnYs;returnZs;returnj;returnk;returnl;}publicvoidMain(string[]args)//计算流程n{while(count=4){calX();calJSZ();calXSJZandCSX();Console.WriteLine(第{0}次迭代的结果,count);calJSGZS();count++;}9}privatevoidMA(double[,]c){inti,j,h,m;constintn=6;doublel;double[,]q=newdouble[n,12];for(i=0;in;i++)//求高斯矩阵{for(j=0;j12;j++){q[i,j]=c[i,j];}}for(i=0;in;i++)//{for(j=0;j12;j++)//构造单位阵{if(i+6==j)q[i,j]=1;elseq[i,j]=0;}}for(h=0,m=0;mn-1;m++,h++)//消去对角线以下的数据{for(i=m+1;in;i++){if(q[i,