摄影测量实习报告

1、.摄影测量实习报告 单张影像空间后方交会程序设计 实习时间 2013.5.20-2013.5.24 学生班级 测绘10-2班 学生姓名 刘 航 学生学号 1072143212 所在院系 矿业工程学院 指导教师 张会战 邵亚琴 1、 实习目的1 深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。2 掌握空间后方交会的定义和实现算法(1) 定义：空间后方交会是以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，。(2) 算法：由

2、于每一对像方和物方共轭点可列出2个方程，因此若有3个已知地面坐标的控制点，则可列出6个方程，解求6个外方位元素的改正数Xs，Ys，Zs，。实际应用中为了提高解算精度，常有多余观测方程，通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点，因而要用最小二乘平差方法进行计算。3. 利用Visual C+或者Matlab（或其他熟悉的计算机语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并进行评定精度。4通过编写程序实现单张影像空间后方交会计算，掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点，巩固各类基础课程及计算机课程的学习

3、内容，培养上机调试程序的动手能力，通过对实验结果的分析，增强综合运用所学知识解决专业实际问题的能力。二、实习环境1硬件环境：Window操作系统2软件环境：VC+或Matlab或其他计算机语言三、实习内容利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素并进行精度评定。4、 实习原理以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素 Xs，Ys，Zs，。共线条件方程如下：其中：x,y为像点的像平面坐标；x0,y0,f为影像的外方位元素；Xs,Ys,Zs为摄站点的物方空间坐标；X,Y,

4、Z为物方点的物方空间坐标；2精度评定 其中输入原始数据归算像点坐标x,y计算和确定初值Xs0, Ys0, Zs0, 0,0,0组成旋转矩阵R计算（x），（y）和lx，ly逐点组成误差方程并法化所有点完否？解法方程，求未知数改正数计算改正后外方位元素未知数改正数<限差否？整理并输出计算结果正常输出 迭代次数小于限差否？输出中间结果和出错信息非正常结束程序框图五、实习数据1点位略图3 9 5 81 76 12 5 84 102模拟像片一对：左片号233像片比例尺： 1/300004航摄机主距：f=150mm5每张像片有4个控制点 6各片像点坐标及其地面坐标片号点号 像点坐标（mm） 地面坐标

5、（m） x y X Y Z23 1 -91.596 -74.859100000.00137500.0011.00 3 -94.230 81.446100000.00142500.0036.00 7 95.207 -75.521106000.00137500.0042.00 9 96.797 83.077106000.00142500.0056.0024 4-102.695 -79.618103000.00137500.0090.00 6 -99.904 81.754103000.00142500.0031.0010 86.890 -77.540109000.00137500.00 7.0012

6、 88.904 76.257109000.00142500.00 5.00六、验证数据1已知航摄仪内方位元素f153.24mm，XoYo0。2已知4对点的影像坐标和地面坐标 影像坐标地面坐标x（mm）y（mm）X(m)Y（m）Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31 3解算结果 4.数据分析 选取第六张像片进行计算，迭代次数为2次。经过比较发现，计算出

7、的6个外方位元素与所给参考值相比，相差很小，计算结果符合要求：线元素误差小于0.5米；角元素误差30秒。计算其精度，可以通过法方程式中未知数的系数矩阵的逆阵A-1来求解，此时，视像点坐标为等精度不相关观测值。线元素精度mx等，高于0.05米；角元素精度高于0.00003弧度。计算结果都达到标准。在此次计算中，我运用了所给的全部控制点，而空间后方交会所运用的控制点，应该避免位于一个圆柱面上，否则会出现解不唯一的情况。选点时，还需要避免选择的点过于聚集在一起，或位于一条直线上，所选控制点最好分布在像片的四角和中央。并且数量充足，这样有利于提高解算精度。迭代时，所选择控制条件不同，迭代次数略有不同，

8、所以最后结果也会略有不同。一般设置为线元素改正数小于0.01m，角元素改正数小于0.1´。 所提供X Y Z为地面测量坐标，带入共线方程时，需要转换为地面摄影测量坐标，最简单的方法为互换XY的数值，即可达到转换坐标目的。并且其单位为米，而像点坐标的单位为厘米，需要统一坐标单位。七、程序输入文件形式如下：C+源程序如下：#include <iostream>#include <fstream>#include <cmath>#include <string>#include <iomanip>using namespace s

9、td;const int n=6;void inverse (double cnn);template<typename T1,typename T2>void transpose (T1*mat1,T2*mat2,int a,int b);template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c);template<typename T>void input (T*mat,int a,int b);template<typen

10、ame T>void output(T*mat,char*s,int a,int b);int main()ofstream outFile;cout.precision(5);double x0=0.0, y0=0.0; double fk=0.15324; /内方位元素double m=39689; /估算比例尺 double B45=0.0,R33,XG61,AT68,ATA66,ATL61;input (B,4,5); /从文件中读取控制点的影像坐标和地面坐标，存入数组Bdouble Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0; double

11、 X,Y,Z,L81,A86; /确定未知数的出始值for(int i=0;i<4;i+)Xs=Xs+Bi2; Ys=Ys+Bi3; Zs=Zs+Bi4;Xs=Xs/4; Ys=Ys/4; Zs=Zs/4+m*fk;int f=0;do/迭代计算f+;/组成旋转矩阵 R00=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K); R01=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K); R02=-sin(Q)*cos(W); R10=cos(W)*sin(K); R11=cos(W)*cos(K); R12=-sin(W); R20=sin(Q)*

12、cos(K)+cos(Q)*sin(W)*sin(K); R21=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K); R22=cos(Q)*cos(W); /计算系数阵和常数项 for(int i=0,k=0,j=0;i<=3;i+,k+,j+) X=R00*(Bi2-Xs)+R10*(Bi3-Ys)+R20*(Bi4-Zs);Y=R01*(Bi2-Xs)+R11*(Bi3-Ys)+R21*(Bi4-Zs);Z=R02*(Bi2-Xs)+R12*(Bi3-Ys)+R22*(Bi4-Zs); Lj0=Bi0-(x0-fk*X/Z);Lj+10=Bi1-(y0-fk*Y

13、/Z);j+;Ak0=(R00*fk+R02*(Bi0-x0)/Z;Ak1=(R10*fk+R12*(Bi0-x0)/Z;Ak2=(R20*fk+R22*(Bi0-x0)/Z;Ak3=(Bi1-y0)*sin(W)-(Bi0-x0)*(Bi0-x0)*cos(K)-(Bi1-y0)*sin(K)/fk+fk*cos(K)*cos(W);Ak4=-fk*sin(K)-(Bi0-x0)*(Bi0-x0)*sin(K)+(Bi1-y0)*cos(K)/fk;Ak5=Bi1-y0;Ak+10=(R01*fk+R02*(Bi1-y0)/Z;Ak+11=(R11*fk+R12*(Bi1-y0)/Z;Ak+

14、12=(R21*fk+R22*(Bi1-y0)/Z;Ak+13=-(Bi0-x0)*sin(W)-(Bi1-y0)*(Bi0-x0)*cos(K)-(Bi1-y0)*sin(K)/fk-fk*sin(K)*cos(W);Ak+14=-fk*cos(K)-(Bi1-y0)*(Bi0-x0)*sin(K)+(Bi1-y0)*cos(K)/fk;Ak+15=-(Bi0-x0);k+;transpose(A,AT,6,8);multi(AT,A,ATA,6,8,6);inverse(ATA);multi(AT,L,ATL,6,8,1);multi(ATA,ATL,XG,6,6,1);Xs=Xs+XG0

15、0; Ys=Ys+XG10; Zs=Zs+XG20;Q=Q+XG30; W=W+XG40; K=K+XG50;while(XG30>=6.0/206265.0|XG40>=6.0/206265.0|XG50>=6.0/206265.0);cout<<"迭代次数为："<<f<<endl;/精度评定double AXG81,V81,VT18,VTV11,m0,D66;multi(A,XG,AXG,8,6,1); for( i=0;i<8;i+) /计算改正数Vi0=AXGi0-Li0; transpose (V,VT,

16、1,8); multi(VT,V,VTV,1,8,1); m0=VTV00/2;for(i=0;i<6;i+)for(int j=0;j<6;j+)Dij=m0*ATAij;/屏幕输出误差方程系数阵、常数项、改正数output(A,"误差方程系数阵A为：",8,6);output(L,"常数项L为：",8,1);output(XG,"改正数为：",6,1);outFile.open("aim.txt",ios:app); /打开并添加aim.txt文件 outFile.precision(10);/以文

17、件的形式输出像片外方位元素、旋转矩阵、方差阵outFile<<"一、像片的外方位元素为："<<endl<<endl;outFile<<setw(10)<<"Xs="<<Xs<<setw(10)<<"Ys="<<Ys<<setw(10)<<"Zs="<<Zs<<endl;outFile<<setw(20)<<"航向倾角为：&quo

18、t;<<Q<<setw(10)<<"旁向倾角为："<<W<<setw(10)<<"像片旋角为："<<K<<endl;outFile<<endl<<"二、旋转矩阵R为:"<<endl<<endl;for( i=0;i<3;i+)for(int j=0;j<3;j+)outFile<<setw(25)<<Rij<<setw(25);outFile&l

19、t;<endl;outFile<<endl;outFile<<setw(0)<<"三、精度评定结果为："<<endl;outFile.precision(5);for(i=0;i<6;i+)for(int j=0;j<6;j+)outFile<<setw(14)<<Dij<<setw(14);outFile<<endl;outFile.close();return 0;template<typename T1,typename T2>void tra

20、nspose(T1*mat1,T2*mat2,int a,int b) int i,j; for(i=0;i<b;i+) for(j=0;j<a;j+)mat2ji=mat1ij;return;template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j,k;for(i=0;i<a;i+)for(j=0;j<c;j+)resultij=0;for(k=0;k<b;k+)resultij+=mat1ik*mat2kj

21、;return;template <typename T>void input (T*mat,int a,int b) ifstream inFile;inFile.open("控制点坐标.txt");while(!inFile.eof()for (int i=0;i<a;i+) for(int j=0;j<b;j+)inFile>>matij;inFile.close();return;template<typename T>void output(T*mat,char*s,int a,int b) cout<<

22、setw(15)<<s<<endl;for(int i=0;i<a;i+)for(int j=0;j<b;j+) cout<<setw(13)<<matij;cout<<endl;return;void inverse(double cnn) int i,j,h,k; double p; double qn12; for(i=0;i<n;i+)/构造高斯矩阵 for(j=0;j<n;j+) qij=cij; for(i=0;i<n;i+) for(j=n;j<12;j+) if(i+6=j) qij=1; else qij=0; for(h=k=0;k<n-1;k+,h+)/消去对角线以下的数据 for(i=k+1;i<n;i+) if(q