單像空間后方交會(huì)實(shí)驗(yàn)報(bào)告（c版)

1、單像空間后方交會(huì)姓名： 學(xué)號(hào)： 時(shí)間： Echo did this for you .2013/4/25目錄一、作業(yè)任務(wù)- 3 -二、計(jì)算原理- 3 -三、算法流程- 7 -四、源程序- 8 -五、計(jì)算結(jié)果- 8 -六、結(jié)果分析- 8 -七、心得與體會(huì)- 8 -八、附頁- 8 -1.c+程序- 8 -2.C+程序截圖- 15 -3.matlb程序- 16 -一、 作業(yè)任務(wù)已知條件：攝影機(jī)主距f=153.24mm，x0=0，y0=0, 像片比例尺為1:40000，有四對(duì)點(diǎn)的像點(diǎn)坐標(biāo)與相應(yīng)的地面坐標(biāo)如下表。 點(diǎn)號(hào)像點(diǎn)坐標(biāo)地面坐標(biāo)x(mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.99

2、36589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31以單像空間后方交會(huì)方法，求解該像片的外方位元素。二、 計(jì)算原理1. 獲取已知數(shù)據(jù)。從航攝資料中查取平均航高與攝影機(jī)主距；獲取控制點(diǎn)的地面測(cè)量坐標(biāo)并轉(zhuǎn)換為地面攝測(cè)坐標(biāo)。2. 測(cè)量控制點(diǎn)的像點(diǎn)坐標(biāo)并作系統(tǒng)誤差改正。3. 確定未知數(shù)的初始值。在豎直攝影且地面控制點(diǎn)大體對(duì)稱分布的情況下，按如下方法確定初始值，即式中：為攝影比例尺分母；為控制點(diǎn)個(gè)數(shù) 4.

3、用三個(gè)角元素的初始值按下式計(jì)算各方向余弦值，組成旋轉(zhuǎn)矩陣矩陣中各元素的計(jì)算公式如下：5. 逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值。利用未知數(shù)的近似值和控制點(diǎn)的地面坐標(biāo)，帶入以下共線方程式， 逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值、6. 逐點(diǎn)計(jì)算誤差方程式的系數(shù)和常數(shù)項(xiàng)，組成誤差方程式。由常數(shù)項(xiàng)計(jì)算公式：得到常數(shù)項(xiàng)矩陣計(jì)算式為：7. 計(jì)算法方程的是系數(shù)矩陣和常數(shù)項(xiàng)，組成法方程式。8. 解法方程，求得外方位元素的改正數(shù)。9. 用前次迭代取得的近似值，加本次迭代的改正數(shù)，計(jì)算外方位元素的新值。式中：代表迭代次數(shù)。10. 將求得的外方位元素改正數(shù)與規(guī)定的限差比較，若小于限差，則迭代結(jié)束。否則用新的近似值重復(fù)49滿足要求為止。11

4、. 由誤差方程的計(jì)算式：以及單位權(quán)中誤差的計(jì)算式（式中：表示多余觀測(cè)數(shù)）以及平差中6個(gè)參數(shù)的協(xié)因數(shù)陣最終得到參數(shù)的中誤差為三、 算法流程否是否是輸入原始數(shù)據(jù)像點(diǎn)坐標(biāo)計(jì)算，系統(tǒng)誤差改正確定外方位元素初始值組成旋轉(zhuǎn)矩陣R逐點(diǎn)組成誤差方程式并法化所有像點(diǎn)完否計(jì)算改正后的外方位元素計(jì)算中誤差，輸出成果，結(jié)束解法方程，求外方位元素改正數(shù)外方位元素改正數(shù)是否小于限差是結(jié)束并提示錯(cuò)誤信息迭代次數(shù)小于n四、 源程序源程序代碼請(qǐng)見附頁五、 計(jì)算結(jié)果在經(jīng)過三次迭代之后得到的最終成果如表1所示表1，最終計(jì)算結(jié)果Xs（米）Ys（米）Zs（米）（弧度）（弧度）（弧度）39795.45225827476.46238575

5、72.685988-0.0039870.0021140.067578六、 結(jié)果分析由計(jì)算結(jié)果可知在拍攝照片瞬間，攝影中心在地面攝影測(cè)量坐標(biāo)系中的坐標(biāo)為（39795.452258，27476.462385，7572.685988）（單位：米），航向傾角為-0.003987弧度，旁向傾角為0.002114弧度，像片旋角為0.067578弧度。表2，精度評(píng)定結(jié)果參數(shù)XsYsZs中誤差1.1073 1.2494 0.4881 0.0002 0.0002 0.0001 七、 心得與體會(huì)通過本次實(shí)驗(yàn)，我對(duì)單張像片的空間后方交會(huì)的計(jì)算原理及實(shí)現(xiàn)過程有了很深的認(rèn)識(shí)，并在牢記各種計(jì)算公式的推導(dǎo)過程基礎(chǔ)上做到了熟

6、練應(yīng)用。在本次實(shí)驗(yàn)中，我遇到了很多困難，但是都一一得到了解決，在解決困難的過程中的編程能力得到了提高，對(duì)其所涉及到的知識(shí)的印象也得到了加深。這個(gè)程序本身還有一些不足，例如:在迭代過程中的矩陣的相乘、輸出等，沒有采用編寫函數(shù)的形式，而是逐一進(jìn)行計(jì)算，這是程序整體看起來較臃腫。在編程過程中為了確保計(jì)算的準(zhǔn)確性，我還利用matlab同步編寫了一個(gè)相同功能的計(jì)算單張相片空間后方交會(huì)的迭代程序，并以此來驗(yàn)算原c+程序的正確性?，F(xiàn)在編寫完這個(gè)程序，心中很自豪，也有一些遺憾，因?yàn)槲蚁胗胏#來編寫，因?yàn)閏#在暑期實(shí)習(xí)用過之后，時(shí)隔近一年在沒有復(fù)習(xí)，不免把學(xué)會(huì)的知識(shí)又還給了老師，本想通過這次實(shí)驗(yàn)對(duì)c#進(jìn)行復(fù)習(xí)，

7、但是由于時(shí)間原因以及一些個(gè)人因素，沒能達(dá)到自己原本的目標(biāo)?？傊?，感謝老師給了我這次鍛煉自己、提升自己的機(jī)會(huì)。八、 附頁1. c+程序#include"iostream"using namespace std;#include"fstream"#include <iomanip>const int n=6;void inverse(double cnn);template<typename T1,typename T2>void transpose (T1*mat1,T2*mat2,int a,int b);template<

8、typename T1,typename T2>void multi(T1*mat1,T2*mat2,T2 * result,int a,int b,int c);int main()double x42=-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,0.01046,0.06443;double X43=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=

9、0;/num為迭代次數(shù)double X06=0;/設(shè)定未知數(shù)(XS,YS,ZS,三個(gè)角度)初始值double f=0.15324;/攝影機(jī)主距f=153.24mmdouble a=1/40000.0;/像片比例尺為1:40000double R33=0;/初始化旋轉(zhuǎn)矩陣R double approx_x8=0;/用于存放像點(diǎn)估計(jì)值。double A86=0;/定義了一個(gè)系數(shù)陣double AT68=0;/定義了A的轉(zhuǎn)置矩陣double L8=0;/定義常數(shù)項(xiàng)const double pi=3.1415926535897932;double Asum66=0;double result26=0;

10、double result168=0;double sumXYZ3=0;cout.precision(5);cout<<"已知像點(diǎn)坐標(biāo)為：n"for(i=0;i<4;i+)for(j=0;j<2;j+)cout<<fixed;if(j=0)cout<<"x"<<i+1<<"= "<<setw(8)<<xij<<" "elsecout<<"y"<<i+1<<

11、;"= "<<setw(6)<<xij<<endl;cout<<"已知地面四個(gè)點(diǎn)的坐標(biāo)為：n"for(i=0;i<4;i+)for(j=0;j<3;j+)if(j=0)cout<<"X"cout<<i+1;cout<<"= "<<Xij<<" "elseif(j=1)cout<<"Y"cout<<i+1;cout<<&qu

12、ot;= "<<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+=Xij;for(i=0;i<2;i+)X0i=sumXYZi/4;/X0,Y0初始化X0i=1/a*f+sumXYZ2/4.0;/對(duì)Z0進(jìn)行初始化doR00=cos(X03)*cos(X

13、05)-sin(X03)*sin(X04)*sin(X05);R01=-cos(X03)*sin(X05)-sin(X03)*sin(X04)*cos(X05);R02=-sin(X03)*cos(X04);R10=cos(X04)*sin(X05);R11=cos(X04)*cos(X05);R12=-sin(X04);R20=sin(X03)*cos(X05)+cos(X03)*sin(X04)*sin(X05);R21=-sin(X03)*sin(X05)+cos(X03)*sin(X04)*cos(X05);R22=cos(X03)*cos(X04);/第一個(gè)像點(diǎn)的估計(jì)值，第一個(gè)點(diǎn)的坐

14、標(biāo)存放于X00,X01,X02approx_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)+R12*(X01-X01)+R22*(X02-X02);/第二個(gè)像點(diǎn)的估計(jì)值,第2個(gè)點(diǎn)的坐標(biāo)存放于X10,X11,X12approx_x2=-f*(R00*(X10-X00)+R10*(X11-X01)+R20*(X12-X02)

15、/(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);/第三個(gè)像點(diǎn)的估計(jì)值,第3個(gè)點(diǎn)的坐標(biāo)存放于X20,X21,X22approx_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

16、)+R11*(X21-X01)+R21*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);/第四個(gè)像點(diǎn)的估計(jì)值,第4個(gè)點(diǎn)的坐標(biāo)存放于X30,X31,X32approx_x6=-f*(R00*(X30-X00)+R10*(X31-X01)+R20*(X32-X02)/(R02*(X30-X00)+R12*(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-X0

17、2);for(i=0;i<4;i+)/第i個(gè)像點(diǎn)估計(jì)值放在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_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)

18、;/*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);/*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*

19、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=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)

20、;/*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;/進(jìn)行常數(shù)項(xiàng)的初始化for(i=0;i<4;i+)L2*i=xi0-approx_x2*i;L2*i+1=xi1-approx_x2*i+1;/A的轉(zhuǎn)置矩陣for(i=0;i<8;i+)for(j=0;j<6;j+)ATji=Aij;/實(shí)現(xiàn)A與AT相乘int k=0;for(i=0;i<6;i+)for(j=0;j<6;j+)As

21、umij=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);/實(shí)現(xiàn)矩陣Asum66與AT68的相乘,結(jié)果存放在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;/實(shí)現(xiàn)result168與l8的相乘，得

22、到結(jié)果放在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+)X0i=X0i+result2i;ofstream f7("d:A.txt");f7<< std:fixed;cout<<"進(jìn)行第"<<num<<"次迭代帶得到Xs,Ys,Zs, ，改正數(shù)分別為：n"for(i=0;i<6;i+)co

23、ut<<setw(12)<<result2i;f7<<setw(12)<<result2i;cout<<endl<<endl;f7.close();getchar();while(abs(result23*206265.0)>6|abs(result24*206265.0)>6|abs(result25*206265.0)>6);cout<<"n滿足限差條件結(jié)束循環(huán)，最終結(jié)果為：n"cout<<setw(12)<<"Xs"<

24、<setw(12)<<"Ys"<<setw(12)<<"Zs"<<setw(12)<<""<<setw(12)<<""<<setw(12)<<""<<endl;ofstream f7("d:A.txt");f7<< std:fixed;cout.precision(4);for(i=0;i<6;i+)cout<<setw(

25、12)<<X0i;f7<<setw(16)<<X0i;f7.close();/今double XG61;for(i=0;i<6;i+)XGi0=result2i;double AXG81,V81,VT18,VTV11,m0,D66;multi(A,XG,AXG,8,6,1); for( i=0;i<8;i+) /計(jì)算改正數(shù)Vi0=AXGi0-Li; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV00/2;cout<<endl;ofstream f6("d:what.tx

26、t");/f6<< std:fixed;for(i=0;i<6;i+)for(int j=0;j<6;j+)Dij=m0*Asumij;cout<<setw(10)<<Dij;f6<<setw(15)<<Dij;cout<<endl;f6<<endl;for(i=0;i<6;i+)cout<<sqrt(Dii)<<endl;f6.close();/屏幕輸出誤差方程系數(shù)陣、常數(shù)項(xiàng)、改正數(shù)/getchar();return 0;void inverse(doubl

27、e cnn) int i,j,h,k;double p;double qn12;for(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;elseqij=0;for(h=k=0;k<n-1;k+,h+)/消去對(duì)角線以下的數(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-) /

28、 消去對(duì)角線以上的數(shù)據(jù)for(i=k-1;i>=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+)/將對(duì)角線上數(shù)據(jù)化為1p=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;template<typename T1,typename T2>void transpose(T1*mat1,T2*mat2,int a,int b) int

29、i,j; for(i=0;i<b;i+) for(j=0;j<a;j+)mat2ji=mat1ij;return;template<typename T1,typename T2>void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j,k;for(i=0;i<a;i+)for(j=0;j<c;j+)resultij=0;for(k=0;k<b;k+)resultij+=mat1ik*mat2kj;return;2. C+程序截圖3. matlb程序function XS,Y

30、S,ZS,phi,w,kappa=doitx_original=zeros(8,1);X=zeros(4,3);x_original=-0.08615;-0.06899;-0.05340;0.08221;-0.01478;-0.07663;0.01046;0.06443;X=36589.41,25273.32,2195.17;37631.08,31324.51,728.69;39100.97,24934.98,2386.50,;40426.54,30319.81,757.31;f=0.15324;%/攝影機(jī)主距f=153.24mma=0.000025;%/像片比例尺為1:40000m=4000

31、0;XS=38437YS=27963.155ZS=m*f+1516.9175phi=0w=0kappa=0R=zeros(3,3);num=1;dx=zeros(6,1);while num=1|abs(dx(4,1)*206265)>6|abs(dx(5,1)*206265)>6|abs(dx(6,1)*206265)>6R(1,1)=cos(phi)*cos(kappa)-sin(phi)*sin(w)*sin(kappa);R(1,2)=-cos(phi)*sin(kappa)-sin(phi)*sin(w)*cos(kappa);R(1,3)=-sin(phi)*co

32、s(w);R(2,1)=cos(w)*sin(kappa);R(2,2)=cos(w)*cos(kappa);R(2,3)=-sin(w);R(3,1)=sin(phi)*cos(kappa)+cos(phi)*sin(w)*sin(kappa);R(3,2)=-sin(phi)*sin(kappa)+cos(phi)*sin(w)*cos(kappa);R(3,3)=cos(phi)*cos(w);x_similar=zeros(8,1);x_similar(1,1)=-f*(R(1,1)*(X(1,1)-XS)+R(2,1)*(X(1,2)-YS)+R(3,1)*(X(1,3)-ZS)/(

33、R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS);%x1x_similar(2,1)=-f*(R(1,2)*(X(1,1)-XS)+R(2,2)*(X(1,2)-YS)+R(3,2)*(X(1,3)-ZS)/(R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS);%y1x_similar(3,1)=-f*(R(1,1)*(X(2,1)-XS)+R(2,1)*(X(2,2)-YS)+R(3,1)*(X(2,3)-ZS)/(R(1,3)*(X(2,1)-XS)+R(2,3)*

34、(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS);%x2x_similar(4,1)=-f*(R(1,2)*(X(2,1)-XS)+R(2,2)*(X(2,2)-YS)+R(3,2)*(X(2,3)-ZS)/(R(1,3)*(X(2,1)-XS)+R(2,3)*(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS);%y2x_similar(5,1)=-f*(R(1,1)*(X(3,1)-XS)+R(2,1)*(X(3,2)-YS)+R(3,1)*(X(3,3)-ZS)/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)

35、-ZS);%x3x_similar(6,1)=-f*(R(1,2)*(X(3,1)-XS)+R(2,2)*(X(3,2)-YS)+R(3,2)*(X(3,3)-ZS)/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)-ZS);%y3x_similar(7,1)=-f*(R(1,1)*(X(4,1)-XS)+R(2,1)*(X(4,2)-YS)+R(3,1)*(X(4,3)-ZS)/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS);%x4x_similar(8,1)=-f*(R(1,2)*(X(4,1)-XS)+R(2,2)*(X(4,2)-YS)+R(3,2)*(X(4,3)-ZS)/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS);%y4a=zeros(8,6);for