后方交会求外方位元素

1、精选优质文档-倾情为你奉上西南交通大学摄影测量学课程设计报告单像空间后方交会求解外方位元素 目录：作业任务3计算原理3算法流程6源程序7计算结果12结果分析13心得体会13作业任务：已知条件摄影机主距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.

2、4664.4340426.5430319.81757.31以单像空间后方交会方法，求解该像片的外方位元素。计算原理：如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。那么由共线条件方程知，其中ai，bi，ci是只含三个独立参数，的九个方向余弦。在方程中共有六个未知参数XS,YS,ZS,，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有用新的符号表示各偏导数后为其中（x）、（

3、y）是函数近似值，,是外方位元素近似值的改正数，它们的系数为函数的偏导数。为了便于推导，令=那么有对于系数，其严密算法（以,为例）如下：对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数,，可列出每个点的误差方程式用矩阵表示为那么由 算法流程：Ø 获取已知数据 m, x0 , y0 , f , Xtp, Ytp, Ztp Ø 量测控制点像点坐标 x,y Ø 确定未知数初值 Xs0, Ys0, Z

4、s0, j0, w0, k0 Ø 组成误差方程式并法化Ø 解求外方位元素改正数Ø 检查迭代是否收敛具体的流程程序框图如下：输入原始数据像点坐标计算，系统误差改正确定外方位元素初始值组成旋转矩阵R结束并显示错误信息是迭代次数小于n逐点组成误差方程式并法化所以像点完否否否解法方程，求外方位元素改正数是计算改正后的外方位元素外方位元素改正数是否小于限差否是输出计算成果，计算结束源程序：程序代码： #include<iostream> #include<cmath> using namespace std;const int N=4; const

5、int n=6; /*-矩阵相乘-*/ void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2) int i,j,k;for(i=0;i<i_1;i+)for(j=0;j<j_2;j+)resulti*j_2+j=0.0;for(k=0;k<j_12;k+)resulti*j_2+j+=m1i*j_12+k*m2j+k*j_2;return;/*-矩阵求逆-*/ void inverse(double cnn) int i,j,h,k; double p; double qn12; for

6、(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(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>

7、=0;i-) if(qih=0) continue; p=qkh/qih; for(j=0;j<12;j+) qij*=p; qij-=qkj;for(i=0;i<n;i+)/将对角线上数据化为1 p=1.0/qii; for(j=0;j<12;j+) qij*=p; for(i=0;i<n;i+) /提取逆矩阵 for(j=0;j<n;j+) cij=qij+6;/*-矩阵转置-*/ void transpose(double *m1,double *m2,int m,int n) /矩阵转置 int i,j; for(i=0;i<m;i+) for(j=

8、0;j<n;j+) m2j*m+i=m1i*n+j; return; void main()double Xs,Ys,Zs,q,w,k;double a3,b3,c3;double x0,y0,f;double xN,yN;double XN,YN,ZN;double x1N,y1N;double m;double L2*N;double XX6;double A2*N6;double X0N,Y0N,Z0N,At62*N,result166,result261;int i,n=0;double sum=0,m0;/*-输入点地面坐标-*/X0=36589.41;X1=37631.08;

9、X2=39100.97;X3=40426.54;Y0=25273.32;Y1=31324.51;Y2=24934.98;Y3=30319.81;Z0=2195.17;Z1=728.69;Z2=2386.50;Z3=757.31;/*-输入点像片坐标-*/x0=-86.15;x1=-53.40;x2=-14.78;x3=10.46;y0=-68.99;y1=82.21;y2=-76.63;y3=64.43; cout<<endl;/*-设定外方位元素初始值-*/x0=0;y0=0;f=153.24;m=40000;Xs=0;Ys=0;Zs=f*m/1000;q=0;w=0;k=0;X

10、X3=1;/*-迭代计算-*/while(XX3>6/ | XX4>6/ | XX5>6/)&&n<100)/*-旋转矩阵R-*/a0=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a1=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a2=-sin(q)*cos(w);b0=cos(w)*sin(k);b1=cos(w)*cos(k);b2=-sin(w);c0=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c1=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k

11、);c2=cos(q)*cos(w);/*-像点坐标计算值-*/for(i=0;i<N;i+)X0i=a0*(Xi-Xs)+b0*(Yi-Ys)+c0*(Zi-Zs);Y0i=a1*(Xi-Xs)+b1*(Yi-Ys)+c1*(Zi-Zs); Z0i=a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs);x1i=x0-f*X0i/Z0i;y1i=y0-f*Y0i/Z0i;/*-误差方程中各偏导数的值-*/for(i=0;i<N;i+)A2*i0=(a0*f+a2*(xi-x0)/Z0i;A2*i1=(b0*f+b2*(xi-x0)/Z0i;A2*i2=(c0*f+c2*

12、(xi-x0)/Z0i;A2*i3=(yi-y0)*sin(w)-(xi-x0)*(xi-x0)*cos(k)-yi*sin(k)/f+f*cos(k) *cos(w);A2*i4=-f*sin(k)-(xi-x0)*(xi-x0)*sin(k)+(yi-y0)*cos(k)/f;A2*i5=yi-y0; L2*i=xi-x1i;A1+2*i0=(a1*f+a2*(yi-y0)/Z0i;A1+2*i1=(b1*f+b2*(yi-y0)/Z0i;A1+2*i2=(c1*f+c2*(yi-y0)/Z0i;A1+2*i3=-(xi-x0)*sin(w)-(yi-y0)*(xi-x0)*cos(k)-

13、(yi-y0)*sin(k)/f-f*sin(k) *cos(w);A1+2*i4=-f*cos(k)-(yi-y0)*(xi-x0)*sin(k)+(yi-y0)*cos(k)/f;A1+2*i5=-xi+x0;L1+2*i=yi-y1i;/*-解法方程-*/transpose(&A00,&At00,2*N,6);mult(&At00,&A00,&result100,6,2*N,6);inverse(result1);mult(&At00,L,&result200,6,2*N,1);mult(&result100,&re

14、sult200,&XX0,6,6,1);Xs+=XX0;Ys+=XX1;Zs+=XX2;q+=XX3;w+=XX4;k+=XX5;n+;/*-旋转矩阵R-*/cout<<"迭代次数为:"<<n<<endl;printf("n像片外方位元素的解n");cout<<"航向顷角q:"<<q<<endl;cout<<"旁向倾角w:"<<w<<endl;cout<<"像片旋角k:"<<k<<endl;cout<<"Xs "<<Xs<<endl; cout<<"Ys "<<Ys<<endl; cout<<"Zs "<<Zs<<endl; cout<<endl; cout