摄影测量后方交会

CopyrightXiayun.Xiang1单张相片后方交会CopyrightXiayun.Xiang2目录作业任务························································3解算原理························································3具体过程························································4算法描述及程序流程··············································4计算结果························································7结果分析························································8心得体会及建议··················································8参考文献························································9CopyrightXiayun.Xiang3一，作业任务已知摄影机主距f=153.24mm，四对点的像点坐标与相应地面坐标列入下表：表1-1计算近似垂直摄影情况下后方交会解。二，解算原理【关键词1】中心投影构像方程在摄影测量学中，最重要的方程就是中心投影构像方程（图2-1）。这个方程将地面点在地面摄影测量坐标系中的坐标（物方坐标）和地面点对应像点的像平面坐标联系起来。在解析摄影测量与数字摄影测量中是极其有用的。在以后将要学习到的双像摄影测量光束法、解析测图仪原理及数字影像纠正等都要用到该式。图2-1在上述公式中：x和y分别为以像主点为原点的像点坐标，相应地面点坐标为X,Y,Z，相片主距f以及外方位元素Xs，Ys，Zs，ψ，ω，κ。而在此次作业中，就是已知四个地面控制点的坐标以及其对应的像点坐标，通过间接平差原理来求解此张航片的外方位元素。【关键词2】间接平差在一个平差问题中，当所选的独立参数X的个数等于必要观测值t时，可将每个观测值表达成这t个参数的函数，组成观测方程，然后依据最小二乘原理求解，这种以观测方程为函数模型的平差方法，就是间接平差方法间接平差的函数模型为：随机模型为：平差准则为：VtPV=min【关键词3】单像空间后方交会利用至少三个已知地面控制点的坐标A（Xa，Ya，Za）、B（Xb，Yb，Zb）、Z（Xc，CopyrightXiayun.Xiang4Yc，Zc），与其影像上对应的三个像点的影像坐标a（xa，ya）、b（xb，yb）、c（xc，yc），根据共线方程，反求该像点的外方位元素Xs，Ys，Zs，ψ，ω，κ。这种解算方法是以单张像片为基础，亦称单像空间后方交会。在此次作业中，就是已知四个控制点在地面摄影测量坐标系中的坐标和对应的像点坐标。由此可以列出8个误差方程，存在两个多余观测数，则n=2。故可利用间接平差里，依据最小二乘法则，进行求解。由于共线条件方程是非线性函数模型，为了便于计算，需把非线性函数表达式用泰勒公式展开成现行形式，即“线性化”。而又因为仅取泰勒级数展开式的一次项，未知数的近似值改正是粗略的，所以必须计算采用逐渐趋近法，解求过程需要反复趋近，直至改正值小于限差为止。三，具体过程1，获取已知点数据：从摄影资料中查取像片比例尺1/m，平均航高，内方元素x0，y0，f；从外业测量成果中，获取控制点的地面测量坐标Xt，Yt，Zt，并转换成摄影测量坐标X,Y,Z。2，量测控制点的像点坐标：将控制点标刺在像片上，利用立体坐标量测仪量测控制点的像框坐标，并经像点坐标改正，得到像点坐标x，y。3，确定未知数的初始值：在竖直摄影测量情况下，角元素的初始值为0，及ψ=ω=κ=0；线元素中，Zso=m*f+(Z[0]+Z[1]+Z[2]+Z[3])/4，Xso，Yso的取值可用四个角点上制点坐标的平均值，即：Xso=(X[0]+X[1]+X[2]+X[3])/4;Yso=(Y[0]+Y[1]+Y[2]+Y[3])/4；4，计算旋转矩阵R：利用角元素的近似值计算方向余弦，组成R阵。公式如下：R[0][0]=cos(ψ)*cos(k)-sin(ψ)*sin(w)*sin(k);R[0][1]=-cos(ψ)*sin(k)-sin(ψ)*sin(w)*cos(k);R[0][2]=-sin(ψ)*cos(w);R[1][0]=cos(w)*sin(k);R[1][1]=cos(w)*cos(k);R[1][2]=-sin(w);R[2][0]=sin(ψ)*cos(k)+cos(ψ)*sin(w)*sin(k);R[2][1]=-sin(ψ)*sin(k)+cos(ψ)*sin(w)*cos(k);R[2][2]=cos(ψ)*cos(w);5，逐点计算像点坐标的近似值：利用未知数的近似值按共线方程计算控制点像点坐标的近似值（x）、（y）；6，组成误差方程式：参照教材（5-8）式、（5-9b）式、（5-4）式逐点计算误差方程的系数阵和常数项。7，组成法方程：计算法方程的系数矩阵与常数项。8，解求外方位元素：根据法方程，按间接平差原理解求外方位元素改正值，并与相应的近似值求和，得到外方位元素的新的近似值。9，检查计算是否收敛：将求得的外方位元素的改正值与规定的限差比较，小于限差则计算终止，否则用新的近似值重复第4至第8步骤计算，直至满足要求为止。四，算法描述及程序流程。算法描述（图4-1）：CopyrightXiayun.Xiang5图4-1程序流程：1，读入数据：新建一个txt格式的文本文档（图4-1），在程序中通过读入数据操作将已知数据导入到程序中。图4-2实现这一操作的语句为：ifstreamshuju;shuju.open(“原始坐标数据.txt”);if(shuju.is_open()){while(!shuju.eof()){for(inti=0;i4;i++){CopyrightXiayun.Xiang6shujux[i];shuju;x[i]=(x[i]-x0)/1000;shujuy[i];shuju;y[i]=(y[i]-y0)/1000;shujuX[i];shuju;shujuY[i];shuju;shujuZ[i];shuju;//分列读入数据}}yuanshi.close();}2，给六个外方位元素赋初值，并计算旋转矩阵R。实现此操作的源程序如下：R[0][0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);R[0][1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);R[0][2]=-sin(q)*cos(w);R[1][0]=cos(w)*sin(k);R[1][1]=cos(w)*cos(k);R[1][2]=-sin(w);R[2][0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);R[2][1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);R[2][2]=cos(q)*cos(w);3，计算系数阵A和常数项L。实现此操作程序如下：for(inti=0,J=0,H=0;i=3;i++){XX=R[0][0]*(X[i]-Xs)+R[1][0]*(Y[i]-Ys)+R[2][0]*(Z[i]-Zs);YY=R[0][1]*(X[i]-Xs)+R[1][1]*(Y[i]-Ys)+R[2][1]*(Z[i]-Zs);ZZ=R[0][2]*(X[i]-Xs)+R[1][2]*(Y[i]-Ys)+R[2][2]*(Z[i]-Zs);xx=-f*XX/ZZ;yy=-f*YY/ZZ;l[J][0]=x[i]-xx;l[J+1][0]=y[i]-yy;J=J+2;A[H][0]=(R[0][0]*f+R[0][2]*x[i])/ZZ;A[H][1]=(R[1][0]*f+R[1][2]*x[i])/ZZ;A[H][2]=(R[2][0]*f+R[2][2]*x[i])/ZZ;A[H][3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[H][4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[H][5]=y[i];A[H+1][0]=(R[0][1]*f+R[0][2]*y[i])/ZZ;A[H+1][1]=(R[1][1]*f+R[1][2]*y[i])/ZZ;A[H+1][2]=(R[2][1]*f+R[2][2]*y[i])/ZZ;A[H+1][3]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[H+1][4]=-f*cos(k)-y[i]/f*(x[i]*sin(k)+y[i]*cos(k));A[H+1][5]=-x[i];H=H+2;}4，根据间接平差原理解算出改正数并与限差（在此程序中，限差是三个角度的改正值的CopyrightXiayun.Xiang7绝对值小于6秒）进行比较，若小于限差则计算完成，否则计算新的外方位元素的近似值，生成新的R阵，系数阵和L阵，重新平差。为实现这一迭代计算，程序中采用了一个dowhile循环语句。While的循环条件是fabs(XF[3][0])*t=6.0||fabs(XF[4][0])*t=6.0||fabs(XF[5][0])*t=6.0。（其中XF[3][0],XF[4][0],XF[5][0]分别为外方位元素中三个角元素的改正值，fabs是取绝对值）2，精度估计及生成解析结果文档（.txt）。另：矩阵相乘及求逆的子函数见附件。五，计算结果2，计算结果：迭代次数为：3次此张相片的外方位元素为：三个线元素：Xs=39795.5m；Ys=27476.5m；Zs=7572.69m三个角元素：ψ=-14.353″；ω=7.61006″；κ=-243.281″精度评定如下：Xs:±1.10739mm；Ys:±1.24952mm；Zs:±0.48813mmψ:±0.000178625″；ω:±0.000161461″；κ:±7.20381e-005″2，生成结果的文本文档如下图（5-1）所示：CopyrightXiayun.Xiang8图5-1程序调试结果如下图（5-2）：图5-1六、结果分析此次单像后方交会的迭代次数为三次，符合要求。可见此次解算六个外方位元素初始值取得很好。通过对外方位元素的最终结果分析，角度近似垂直，并且线元素改正值也不大，所以此张相片的摄影为近似垂直摄影。通过精度评估发现，此次观测角度的精度远高于线元素的观测精度。七，心得体会及建议1，“尽信书则不如无书”，此次编程中我发现了书上(摄影测量学·武汉大学出版社·2005年9月)的两个错误的公式，一个是P66页上关于Zso确定的公式。书上的是Zso=m*f，而正确的是Zso=m*f+(Z[0]+Z[1]+Z[2]+Z[3])/4。还有一个是P66页关于系数阵的一个公式。书上的公式是a24=-xsinw-【x*(x*cosk-ysink)/f-f*sink】cosw，而正确的公式应该是a24=-xsinw-【y*(x*cosk-ysink)/f-f*sink】cosw。由于这两个公式的错误，给我的编程造成了很大的困难。后来才发现这两个公式有误，程序才能继续进行下去。2，闭门造车是造不出好车的，要集思广益。此次编程过程中，我与测绘一班的李志鹏同学积极讨论，请教学长。后来通过分析张勇学长的程序，借鉴了张勇学长关于系数阵和L阵的构造。当时自己没有想到通过数组来完成这个操作，只