单向后方交会实验报告

1、精选优质文档-倾情为你奉上课程设计报告 班级：测绘一班 学号: 姓名：罗超 日期：2015.4.26目录4333一、计算原理 已知条件摄影机主距f=153.24mm，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。 点号像点坐标地面坐标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.3

2、1以单像空间后方交会方法，求解该像片的外方位元素。二、算法流程（1）获取已知数据。从航摄资料中差取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影坐标。（2）量测控制点的像点坐标并作系统误差改正。（3）确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即，（4） 用三个角元素的初始值按下式，计算各个方向余弦值，组成旋转矩阵R （5） 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标；带入共线方程式，逐点近似像点坐标的近似值（x）、（y）。（6） 逐点计算误差方程式的系数和常数项，组成误差方程式。（7） 计算法方程的系数矩阵和常数项

3、，组成法方程式。 （8） 解法方程，求得外方位元素的改正数。（9） 用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。（10） 将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。负责用新的近似值重复（4）-（9），直到满足要求为止。 用共线方程进行空间后方交会的程序框如图所示。输入原始数据像点坐标计算，系统误差正确定外方位因素初始值组成旋转矩阵R逐点组成误差方程式并法化所有像点完否是解法方程，求外方位元素改正数计算改正后的外方位元素外方位元素改正数是否小于限差是输出计算成果，计算并结束结束并显示错误信息迭代次数小于n 是否 否3、 源程序#include”ios

4、tream"using namespace std;#include”fstream"#include <iomanip>const int n=6;void inverse(double cnn);template<typename Tl,typename T2>void transpose (TI*mat1,T2*mat2,int a,int b);template<typename TI,typename T2>void multi(T1*mat1,T2*mat2,T2*result, int a, int b, int c);in

5、t main() doublex42=-0.08615,-0.06899,-0.05340,0.08221，-0.01478，-0.07663,0.01046,0.06443; doubleX43=36589.41,25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98, 2386.50,40426.54,30319.81,757.31; int i,j,num=0;/num为迭代次数 double X06=0;/设定未知数(XS,YS,ZS,三个角度)初始值 double f=0.15324;/摄影机主距f=153.24mm d

6、ouble a=1/40000.0;/像片比例尺为1:40000 double R33=0;/初始化旋转知阵R double approx_ x8=0;/用于存放像点估计值 double A86二0;/定义了一个系数阵 double AT68二0;/定义了A的转置知阵 double L8=0:/定义常数项const double pi=3.97932;double Asum66=0;double result26=0;double resultl68=0;double sumXYZ3=0;cout.precision(5);cout<<”已知像点坐标为:n”;for(i=0;i&l

7、t;4;i+) for(j=0;j<2;j+) cout<<fixed: if (j=0) cout<<”x”<<i+l<<”=”<<setw (8)<<xij<<” ”;elsecout<<”y"<<i+1<<”=”<<setw(6)<<xij<<endl; cout<<”己知地面四个点的坐标为:n”;for(i=0;i<4;i+) for(j=0;j<3;j+) if (j=0)cout<&

8、lt;”X"cout<<i+1;cout<<”= “<<Xij<<” ”;else if(j=1) cout<<”Y"cout<<i+1;cout<<”= “<<Xij)<<” ”:elsecout<<”Z”;cout<<i+1;cout<<”= ”;cout<<Xij<<endl; cout<<endl;for(j=0;j<3;j+) for(i=0;i<4;i+) sumXYZj+

9、=Xij;for(i=0;i<2;i+) X0i=sumXYZi/4;/X0,Y0初始化X0i=1/a*f+sumXYZ2/4. 0;/对Z0进行初始化do R00=cos(X03)*cos(X05)-sin(X03)*sin(X04)*sin(X05); R0I=-cos(X03)*sin(X05)-sin(X03)*sin(X04)*cos(X05); R02=-sin(X03)*cos(X04); R10=cos(X04)*sin(X05); RI11=cos(X04)*cos(X05); R12=-sin(X04); R20=sin(X03)*cos(X05)+cos(X03)*

10、sin(X04)*sin(X05); R2l=-sin(X03)*sin(X05)+cos(X03)*sin(X04)*cos(X05); R22=cos(X03)*cos(X04);/第一个像点的估计值，第一个点的坐标存放于X 0 0,X 0 1,X 0 2approx_x0=-f*(R00*(X00-X00)+R10*(X01-X01)+R20*(X02-X02)/(R02*(X00-X00)+R12*(X01-X01)+R22*(X02-X02); approx_x1=-f*(R01*(X00-X00)+R11*(X01-X01)+R21*(X02-X02)/(R02*(X00-X00)

11、+R12*(X01-X01)+R22*(X02-X02); /第二个像点的估计值,第2个点的坐标存放于X10,X11,X12 approx_x2=-f*(R00*(X10-X00)+R10*(X11-X01)+R20*(X12-X02)/(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02); approx_x3=-f*(R01*(X10-X00)+R11*(X11-X01)+R21*(X12-X02)/(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02); /第三个像点的估计值,第3个点的坐标存放于X20,X21,X22 appr

12、ox_x4=-f*(R00*(X20-X00)+R10*(X21-X01)+R20*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);approx_x5=-f*(R01*(X20-X00)+R11*(X21-X01)+R21*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02); /第四个像点的估计值,第4个点的坐标存放于X30,X31,X32 approx_x6=-f*(R00*(X30-X00)+R10*(X31-X01)+R20*(X32-X02)/(R02*(X30-X00)+R

13、12*(X31-X01)+R22*(X32-X02); approx_x7=-f*(R01*(X30-X00)+R11*(X31-X01)+R21*(X32-X02)/(R02*(X30-X00)+R12*(X31-X01)+R22*(X32-X02); for(i=0;i<4;i+) /第i个像点估计值放在approx_x2*(i-1),approx_x2*i-1 /*a11*/A2*i0=(R00*f+R02*approx_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02); /*a12*/A2*i1=(R10*f+R12*approx_

14、x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02); /*a13*/A2*i2=(R20*f+R22*approx_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02); /*a21*/A2*i+10=(R01*f+R02*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02); /*a22*/A2*i+11=(R11*f+R12*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);

15、 /*a23*/A2*i+12=(R21*f+R22*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02); /*a14*/A2*i3=approx_x2*i+1*sin(X04)-(approx_x2*i/f*(approx_x2*i*cos(X05)-approx_x2*i+1*sin(X05)+f*cos(X05)*cos(X04); /*a15*/A2*i4=-f*sin(X05)-approx_x2*i/f*(approx_x2*i*sin(X05)+approx_x2*i+1*cos(X05); /*a16*/A2*i5

16、=approx_x2*i+1; /*a24*/A2*i+13=-1*approx_x2*i*sin(X04)-(approx_x2*i+1/f*(approx_x2*i*cos(X05)-approx_x2*i+1*sin(X05)-f*sin(X05)*cos(X04); /*a25*/A2*i+14=-1*f*cos(X05)-approx_x2*i+1/f*(approx_x2*i*sin(X05)+approx_x2*i+1*cos(X05); /*a26*/A2*i+15=-approx_x2*i; /进行常数项的初始化for(i=0;i<4;i+) L2*i=xi0-appr

17、ox_x2*i; L2*i+1=xi1-approx_x2*i+1; /A的转置矩阵 for(i=0;i<8;i+) for(j=0;j<6;j+) ATji=Aij; /实现A与AT相乘 int k=0; for(i=0;i<6;i+) for(j=0;j<6;j+) Asumij=0; for(i=0;i<6;i+) for(k=0;k<6;k+) for(j=0;j<8;j+) Asumik+=ATij*Ajk; /得到AT*A的逆矩阵存放在inverseAsum66中 inverse(Asum); /实现矩阵Asum66与AT68的相乘,结果存

18、放在result168中 for(i=0;i<6;i+) for(j=0;j<8;j+) result1ij=0; for(i=0;i<6;i+) for(k=0;k<8;k+) for(j=0;j<6;j+) result1ik+=Asumij*ATjk; /实现result168与l8的相乘，得到结果放在result26中； for(i=0;i<6;i+) result2i=0; for(i=0;i<6;i+) for(j=0;j<8;j+) result2i+=result1ij*Lj; num+; for(i=0;i<6;i+) X

19、0i=X0i+result2i; ofstream f7("d:A.txt"); f7<< std:fixed; cout<<"进行第"<<num<<"次迭代带得到Xs,Ys,Zs, ，改正数分别为：n" for(i=0;i<6;i+) cout<<setw(12)<<result2i;f7<<setw(12)<<result2i; cout<<endl<<endl; f7.close(); getchar()

20、; while(abs(result23*.0)>6|abs(result24*.0)>6|abs(result25*.0)>6); cout<<"n满足限差条件结束循环，最终结果为：n" cout<<setw(12)<<"Xs"<<setw(12)<<"Ys"<<setw(12)<<"Zs"<<setw(12)<<""<<setw(12)<<&q

21、uot;"<<setw(12)<<""<<endl; ofstream f7("d:A.txt"); f7<< std:fixed; cout.precision(4); for(i=0;i<6;i+) cout<<setw(12)<<X0i;f7<<setw(16)<<X0i; f7.close(); /今 double XG61; for(i=0;i<6;i+) XGi0=result2i; double AXG81,V81,VT18

22、,VTV11,m0,D66; multi(A,XG,AXG,8,6,1); for( i=0;i<8;i+) /计算改正数 Vi0=AXGi0-Li; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV00/2; cout<<endl; ofstream f6("d:what.txt"); / f6<< std:fixed; for(i=0;i<6;i+) for(int j=0;j<6;j+) Dij=m0*Asumij; cout<<setw(10)<<

23、Dij; f6<<setw(15)<<Dij; cout<<endl; f6<<endl; for(i=0;i<6;i+) cout<<sqrt(Dii)<<endl;f6.close(); /屏幕输出误差方程系数阵、常数项、改正数 / getchar(); return 0; 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=

24、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(qih=0) continue; p=qkh/qih; for(j=0;j<12;j+) qij*=p; qij-=qkj; for(h=k=n-1;k>0;k-,h-) / 消去对角线以上的数据 for(i=k-1;i>=0;i-) if(qih=0) continue; p=qkh/qih; for(j=0;j<12;j+) qij*=p; qij-=qkj;