后方交會(huì)求外方位元素

1、精選優(yōu)質(zhì)文檔-傾情為你奉上西南交通大學(xué)攝影測量學(xué)課程設(shè)計(jì)報(bào)告單像空間后方交會(huì)求解外方位元素 目錄：作業(yè)任務(wù)3計(jì)算原理3算法流程6源程序7計(jì)算結(jié)果12結(jié)果分析13心得體會(huì)13作業(yè)任務(wù)：已知條件攝影機(jī)主距f=153.24mm，x0=0，y0=0, 像片比例尺為1:40000，有四對點(diǎn)的像點(diǎn)坐標(biāo)與相應(yīng)的地面坐標(biāo)如下表。 點(diǎn)號像點(diǎn)坐標(biāo)地面坐標(biāo)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以單像空間后方交會(huì)方法，求解該像片的外方位元素。計(jì)算原理：如圖所示，物點(diǎn)A和攝影中心S在地面攝影測量坐標(biāo)系中的坐標(biāo)依次是（X,Y,Z）、（XS,YS,ZS）；像點(diǎn)a在像空間坐標(biāo)系中的坐標(biāo)是（x,y,-f）。那么由共線條件方程知，其中ai，bi，ci是只含三個(gè)獨(dú)立參數(shù)，的九個(gè)方向余弦。在方程中共有六個(gè)未知參數(shù)XS,YS,ZS,，所以有三個(gè)不在一條直線上的已知地面點(diǎn)坐標(biāo)就可以求出像片的這六個(gè)外方位元素。由于共線條件方程是非線性方程，為了便于迭代計(jì)算，需要把方程用泰勒級數(shù)展開，取一次項(xiàng)得到線型表達(dá)式，有用新的符號表示各偏導(dǎo)數(shù)后為其中（x）、（

3、y）是函數(shù)近似值，,是外方位元素近似值的改正數(shù)，它們的系數(shù)為函數(shù)的偏導(dǎo)數(shù)。為了便于推導(dǎo)，令=那么有對于系數(shù)，其嚴(yán)密算法（以,為例）如下：對于豎直攝影而言，像片的角方位元素都是小值，因而各系數(shù)的近似值為為了提高精度和可靠性，通常需要測量四個(gè)或更多的地面控制點(diǎn)和對應(yīng)的像點(diǎn)坐標(biāo)，采用最小二乘平差方法解算。此時(shí)像點(diǎn)坐標(biāo)（x,y）作為觀測值，加入相應(yīng)的偶然誤差改正數(shù),，可列出每個(gè)點(diǎn)的誤差方程式用矩陣表示為那么由 算法流程：Ø 獲取已知數(shù)據(jù) m, x0 , y0 , f , Xtp, Ytp, Ztp Ø 量測控制點(diǎn)像點(diǎn)坐標(biāo) x,y Ø 確定未知數(shù)初值 Xs0, Ys0, Z

4、s0, j0, w0, k0 Ø 組成誤差方程式并法化Ø 解求外方位元素改正數(shù)Ø 檢查迭代是否收斂具體的流程程序框圖如下：輸入原始數(shù)據(jù)像點(diǎn)坐標(biāo)計(jì)算，系統(tǒng)誤差改正確定外方位元素初始值組成旋轉(zhuǎn)矩陣R結(jié)束并顯示錯(cuò)誤信息是迭代次數(shù)小于n逐點(diǎn)組成誤差方程式并法化所以像點(diǎn)完否否否解法方程，求外方位元素改正數(shù)是計(jì)算改正后的外方位元素外方位元素改正數(shù)是否小于限差否是輸出計(jì)算成果，計(jì)算結(jié)束源程序：程序代碼： #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+)/構(gòu)造高斯矩陣 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+)/消去對角線以下的數(shù)據(jù) 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-) / 消去對角線以上的數(shù)據(jù) 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+)/將對角線上數(shù)據(jù)化為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;/*-矩陣轉(zhuǎn)置-*/ void transpose(double *m1,double *m2,int m,int n) /矩陣轉(zhuǎ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;/*-輸入點(diǎn)地面坐標(biāo)-*/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;/*-輸入點(diǎn)像片坐標(biāo)-*/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;/*-設(shè)定外方位元素初始值-*/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;/*-迭代計(jì)算-*/while(XX3>6/ | XX4>6/ | XX5>6/)&&n<100)/*-旋轉(zhuǎn)矩陣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);/*-像點(diǎn)坐標(biāo)計(jì)算值-*/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;/*-誤差方程中各偏導(dǎo)數(shù)的值-*/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+;/*-旋轉(zhuǎn)矩陣R-*/cout<<"迭代次數(shù)為:"<<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