空間后方交會編程

1、學西南交攝影測量學課程設計報告單片空間后方交匯年級:2009級專業(yè):遙感科學與技術分數(shù):教師簽名二零一一年四目錄TOC o 1-5 h z HYPERLINK l bookmark4作業(yè)任務1. HYPERLINK l bookmark6計算原理1. HYPERLINK l bookmark8算法流程2. HYPERLINK l bookmark10源程序3. HYPERLINK l bookmark12計算結(jié)果5. HYPERLINK l bookmark14結(jié)果分析6.心得體會錯.誤！未定義書簽。西南交通大學本科課程設計（論文）第西南交通大學本科課程設計（論文）第1頁作業(yè)任務已知條件攝影機

2、主距f=153.24mm,xO=O,yO=O,像片比例尺為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.31以單像空間后方交會方法,求解該像片的外方位元素。計算原理因為本算法中涉及眾多的矩陣計算，若用C+,C#之類的語言的話在矩陣計算方面就比較

3、麻煩，而眾所周知的matlab（矩陣實驗室）在處理矩陣計算方面的便捷上有著本質(zhì)上的優(yōu)勢。所以本程序設計用matlab來實現(xiàn)。首先通過讀取.txt文件獲得坐標點的信息，并計算得到初始的六個外方位元素，注意其中的Z應該取四個坐標值的平均再加上m*f。然后通過共線方程計算出（x）和（y）,并進一步求得1矩陣。然后組成系數(shù)矩陣A。利用間接平差公式算的改正數(shù)為Al=inv（A*A）*A*1.檢驗得到的改正數(shù)是否符合精度，若不符合就在用算出的改正數(shù)加上計算前的初始值再進行計算直到精度合格為止西南交通大學課程設計(論文)第西南交通大學課程設計(論文)第2頁西南交通大學課程設計(論文)第西南交通大學課程設計(

4、論文)第2頁算法流程開始否結(jié)束判斷內(nèi)方位元素改正值是否收斂/讀取文件獲得四個點的像點坐標和地面/坐標。并計算得外方位元素的預設值。計算改正后的外方位元素組成矩陣AA(8,6)和1(8,1)計算內(nèi)方位元素的改正值是輸出外方位元素定義旋轉(zhuǎn)矩陣R西南交通大學課程設計（論文）第西南交通大學課程設計（論文）第 頁源程序本程序中我編寫了兩個文件，其中一個為主文件(dushu.m)(類似C+中的main函數(shù))另一個寫的是求解改正數(shù)的函數(shù)(bianhuan.m)。第一個函數(shù)(dushu.m)：%讀入文件shuju.txt并將里面的數(shù)據(jù)存入矩陣A【5.4】中fid=fopen(shuju.txt,r);A=fs

5、canf(fid,%f,5,4);%將像點坐標轉(zhuǎn)換為m為單位的數(shù)據(jù)A(1:2,:)=A(1:2,:)/1000;%計算外方位元素的初始值，組成初始值數(shù)組Xx0=mean(A(3:5,:);x0(3)=x0(3)+0.15324*40000;X=x0,000;%第一次運算x0=bianhuan(A,X);%檢驗增量是否合格p=0;i=1;while(p=all(x0(4:6,1)*206265-6)X=x0+X;x0=bianhuan(A,X);i=i+1enddisp(循環(huán)次數(shù))disp(i)disp(外方位元素精度)fori=1:6disp(x0(i,1)enddisp(外方位元素)fori

6、=1:6disp(X(i)end第二個函數(shù)(bianhuanm):functionAl=bianhuan(A,AO%A為原始度數(shù)中的數(shù)組A0為估計后的數(shù)組%定義fai=sin，fai1=cos，kama=sink,kama1=cosk,omig=sinwomig1=coswfai=sin(A0(4);fai1=cos(A0(4);kama=sin(A0(6);kama1=cos(A0(6);omig=sin(A0(5);omig1=cos(A0(5);%定義R矩陣中的a1a2a3b1b2b3c1c2c3a1=fai1*kama1-fai*omig*kama;a2=-fai1*kama-fai*

7、omig*kama1;a3=-fai*omig1;b1=omig1*kama;b2=omig1*kama1;b3=-omig;c1=fai*kama1+fai1*omig*kama;c2=-fai*kama+fai1*omig*kama1;c3=fai1*omig1;%寫出R1(原本的R)矩陣并添設一個R矩陣供后面的計算用R=a1b1c1;a2b2c2;a3b3c3;R1=a1a2a3;b1b2b3;c1c2c3;%將四個坐標的(X-Xs)(Y-Ys)(Z-Zs)設成一個4*3的矩陣方便計算z=zeros(4,3);z(1,:)=(A(3,1)-A0(1)(A(4,1)-A0(2)(A(5,1

8、)-A0(3);z(2,:)=(A(3,2)-A0(1)(A(4,2)-A0(2)(A(5,2)-A0(3);z(3,:)=(A(3,3)-A0(1)(A(4,3)-A0(2)(A(5,3)-A0(3);z(4,:)=(A(3,4)-A0(1)(A(4,4)-A0(2)(A(5,4)-A0(3);%計算1。因為l=【lx%ly】所以用i控制坐標點，用j控制l數(shù)組的位置。其中因為lx=x+f*a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)/a3*(X-X%s)+b3*(Y-Ys)+c3*(Z-Zs),其中的al*(X-Xs)+bl*(Y-Ys)+cl*(Z-Zs)/a3*(X-X%s

9、)+b3*(Y-Ys)+c3*(Z-Zs)換位了矩陣點乘并求和(因為他們都是一維數(shù)組)sum(R(l,:).*z(i,:)/sum(R(3,:).*z(i,:)j=l;fori=l:4l(j,l)=A(l,i)+0.l5324*sum(R(l,:).*z(i,:)/sum(R(3,:).*z(i,:);j=j+l;l(j,l)=A(2,i)+0.l5324*sum(R(2,:).*z(i,:)/sum(R(3,:).*z(i,:);j=j+l;end%定義系數(shù)矩陣AA。其中i控制坐標點，j控制系數(shù)矩陣中的位置AA=zeros(8,6);forj=2:2:8i=j/2;all=l/sum(R(3

10、,:).*z(i,:)*(al*0.l5324+a3*A(l,i);al2=l/sum(R(3,:).*z(i,:)*(bl*0.l5324+b3*A(l,i);al3=l/sum(R(3,:).*z(i,:)*(cl*0.l5324+c3*A(l,i);al4=A(2,i)*omig-(A(l,i)/0.l5324*(A(l,i)*kamal-A(2,i)*kama)+0.l5324*kamal)*omigl;a24=-A(l,i)*omig-(A(2,i)/0.l5324*(A(l,i)*kamal-A(2,i)*kama)-0.l5324*kama)*omigl;al5=-0.l5324

11、*kama-A(l,i)/0.l5324*(A(l,i)*kama+A(2,i)*kamal);al6=A(2,i);a21=1/sum(R(3,:).*z(i,:)*(a2*0.15324+a3*A(2,i);a22=1/sum(R(3,:).*z(i,:)*(b2*0.15324+b3*A(2,i);a23=1/sum(R(3,:).*z(i,:)*(c2*0.15324+c3*A(2,i);a25=-0.15324*kama1-A(2,i)/0.15324*(A(1,i)*kama+A(2,i)*kama1);a26=-A(1,i);j=j-1;AA(j,:)=a11a12a13a14a15a1