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