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;