單像空間后方交會程序《Java版》

1、本科實驗報告 課程名稱：單像空間后方交會編程實現(xiàn) 姓 名：李恒松 學(xué) 院：地學(xué)學(xué)院 專 業(yè)：測繪工程一班 學(xué) 號：20133294 指導(dǎo)老師：陳強 2015年 5 月 20 日目錄一、作業(yè)任務(wù)- 3 -二、計算原理- 3 -三、算法流程- 7 -四、源程序- 8 -五、計算結(jié)果- 8 -六、結(jié)果分析- 8 -七、心得與體會- 8 -八、附頁- 8 -Java程序- 8 -一、 作業(yè)任務(wù)已知條件：攝影機主距f=153.24mm，x0=0，y0=0, 像片比例尺為1:40000，有四對點的像點坐標(biāo)與相應(yīng)的地面坐標(biāo)如下表。 點號像點坐標(biāo)地面坐標(biāo)x(mm)y(mm)X(m)Y(m)Z(m)1-86.1

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

3、點個數(shù) 4. 用三個角元素的初始值按下式計算各方向余弦值，組成旋轉(zhuǎn)矩陣其中：5. 逐點計算像點坐標(biāo)的近似值。利用未知數(shù)的近似值和控制點的地面坐標(biāo)，帶入以下共線方程式， 逐點計算像點坐標(biāo)的近似值、，將題目所給數(shù)據(jù)代入上式可得像點坐標(biāo)的近似數(shù)據(jù)。6. 逐點計算誤差方程式的系數(shù)和常數(shù)項，組成誤差方程式。其中其他項類似可得。由常數(shù)項計算公式：將題目所給數(shù)據(jù)代入上式，可得常數(shù)項。7. 計算法方程的是系數(shù)矩陣和常數(shù)項，組成法方程式。8. 解法方程，求得外方位元素的改正數(shù)。9. 用前次迭代取得的近似值，加本次迭代的改正數(shù)，計算外方位元素的新值。式中：代表迭代次數(shù)。10. 將求得的外方位元素改正數(shù)與規(guī)定的限差

4、比較，若小于限差，則迭代結(jié)束。否則用新的近似值重復(fù)49滿足要求為止。11. 由誤差方程的計算式：以及單位權(quán)中誤差的計算式（式中：表示多余觀測數(shù)）以及平差中6個參數(shù)的協(xié)因數(shù)陣最終得到參數(shù)的中誤差為一、 算法流程否是否是輸入原始數(shù)據(jù)像點坐標(biāo)計算，系統(tǒng)誤差改正確定外方位元素初始值組成旋轉(zhuǎn)矩陣R逐點組成誤差方程式并法化所有像點完否計算改正后的外方位元素計算中誤差，輸出成果，結(jié)束解法方程，求外方位元素改正數(shù)外方位元素改正數(shù)是否小于限差是結(jié)束并提示錯誤信息迭代次數(shù)小于n二、 源程序源程序代碼請見附頁三、 計算結(jié)果在經(jīng)過三次迭代之后得到的最終成果如表1所示表1，最終計算結(jié)果Xs（米） Ys（米）Zs（米）（

5、弧度）（弧度）（弧度）39795.4522575952727476.4623852405647572.685987778016-0.0039869235599752880.00211388806645731-0.06757797053971112四、 結(jié)果分析由計算結(jié)果可知在拍攝照片瞬間，攝影中心在地面攝影測量坐標(biāo)系中的坐標(biāo)為（39795.452258，27476.462385，7572.685988）（單位：米），航向傾角為-0.003987弧度，旁向傾角為0.002114弧度，像片旋角為0.067578弧度。表2，精度評定結(jié)果參數(shù)中誤差Xs1.15327223063957Ys1.30124

6、4230747318Zs0.50838032661769861.8602199035719645E-41.6814981476239593E-47.501994950982504E-5五、 心得與體會通過這次大作業(yè)，我的心得體會以下幾點：第一：編程確實是一件非誠艱難的事情，即使整個單像空間后方交會實驗的原理非常簡單，清晰明了，但是要用編程把它表示出來還是一件非常艱難的事情。我遇到的最大的難點是在矩陣求逆的步驟上，為了解決這個問題。我用了兩天的時間，還請教了數(shù)學(xué)學(xué)院的同學(xué)，程序里求逆矩陣的方法是來自那位同學(xué)的。這次編程也讓我對于線性代數(shù)中的求逆矩陣的問題理解的更加深刻。第二：在發(fā)現(xiàn)問題，解決問題

7、的過程中，我也認(rèn)識到了自己基礎(chǔ)知識的不足，對于Java語言的不熟悉讓我在很多常見的bug面前無能為力。以后要加強編程語言的學(xué)習(xí)。第三：很多東西看似簡單，其實細(xì)節(jié)處才體現(xiàn)出能力。把小的細(xì)節(jié)處理好，才能有大的作品。六、 附頁1. Java程序以下第一部分為主程序import java.util.*;import java.math.*;public class li public static void main(String args)heng hengFind=new heng();double x=-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.

8、07663,0.01046,0.06443;double 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;int i,j,num=0;double X0=new double6;double f=0.15324;double a=1/40000.0;double R=new double33;double evaluate_x=new double8;double A=new double86;double A_T=new double

9、68;double L=new double8;double pi=3.1415926535897932;double AT_A=new double66;double inverseAT_A=new double66;/inverseAT_A=hengFind.inverse(AT_A,6);double result_v=new double8;/用來存放用于計算單位權(quán)中誤差的v值double standard=new double6;double resultF=0;/用來存放v的轉(zhuǎn)置矩陣double result1=new double68;double result2=new dou

10、ble6;double sumXYZ=new double3;System.out.println("已知像點坐標(biāo)為： ");for(i=0;i<4;i+)for(j=0;j<2;j+)if(j=0)System.out.println("x"+(i+1)+"= "+"t"+xij);elseSystem.out.println("y"+(i+1)+"= "+"t"+xij);System.out.println("已知地面四個點的坐

11、標(biāo)為：");for(i=0;i<4;i+)for(j=0;j<3;j+)if(j=0)System.out.println("X"+(i+1)+"= "+"t"+Xij);elseif(j=1)System.out.println("Y"+(i+1)+"= "+"t"+Xij);elseSystem.out.println("Z"+(i+1)+"= "+"t"+Xij);for(j=0;j<

12、;3;j+)for(i=0;i<4;i+)sumXYZj+=Xij;for(i=0;i<2;i+)X0i=sumXYZi/4;X0i=1/a*f+sumXYZ2/4.0;doR00=Math.cos(X03)*Math.cos(X05)-Math.sin(X03)*Math.sin(X04)*Math.sin(X05);R01=-Math.cos(X03)*Math.sin(X05)-Math.sin(X03)*Math.sin(X04)*Math.cos(X05);R02=-Math.sin(X03)*Math.cos(X04);R10=Math.cos(X04)*Math.si

13、n(X05);R11=Math.cos(X04)*Math.cos(X05);R12=-Math.sin(X04);R20=Math.sin(X03)*Math.cos(X05)+Math.cos(X03)*Math.sin(X04)*Math.sin(X05);R21=-Math.sin(X03)*Math.sin(X05)+Math.cos(X03)*Math.sin(X04)*Math.cos(X05);R22=Math.cos(X03)*Math.cos(X04);evaluate_x0=-f*(R00*(X00-X00)+R10*(X01-X01)+R20*(X02-X02)/(R0

14、2*(X00-X00)+R12*(X01-X01)+R22*(X02-X02);evaluate_x1=-f*(R01*(X00-X00)+R11*(X01-X01)+R21*(X02-X02)/(R02*(X00-X00)+R12*(X01-X01)+R22*(X02-X02);evaluate_x2=-f*(R00*(X10-X00)+R10*(X11-X01)+R20*(X12-X02)/(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02);evaluate_x3=-f*(R01*(X10-X00)+R11*(X11-X01)+R21*(X12-X02)/

15、(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02);evaluate_x4=-f*(R00*(X20-X00)+R10*(X21-X01)+R20*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);evaluate_x5=-f*(R01*(X20-X00)+R11*(X21-X01)+R21*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);evaluate_x6=-f*(R00*(X30-X00)+R10*(X31-X01)+R20*(X32-X0

16、2)/(R02*(X30-X00)+R12*(X31-X01)+R22*(X32-X02);evaluate_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+)/*a11*/A2*i0=(R00*f+R02*evaluate_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a12*/A2*i1=(R10*f+R12*evaluate_x2*i)/(R02*(Xi0-X00)

17、+R12*(Xi1-X01)+R22*(Xi2-X02);/*a13*/A2*i2=(R20*f+R22*evaluate_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a21*/A2*i+10=(R01*f+R02*evaluate_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a22*/A2*i+11=(R11*f+R12*evaluate_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a23*/A2*i+12=(R

18、21*f+R22*evaluate_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a14*/A2*i3=evaluate_x2*i+1*Math.sin(X04)-(evaluate_x2*i/f*(evaluate_x2*i*Math.cos(X05)-evaluate_x2*i+1*Math.sin(X05)+f*Math.cos(X05)*Math.cos(X04);/*a15*/A2*i4=-f*Math.sin(X05)-evaluate_x2*i/f*(evaluate_x2*i*Math.sin(X05)+evaluat

19、e_x2*i+1*Math.cos(X05);/*a16*/A2*i5=evaluate_x2*i+1;/*a24*/A2*i+13=-1*evaluate_x2*i*Math.sin(X04)-(evaluate_x2*i+1/f*(evaluate_x2*i*Math.cos(X05)-evaluate_x2*i+1*Math.sin(X05)-f*Math.sin(X05)*Math.cos(X04);/*a25*/A2*i+14=-1*f*Math.cos(X05)-evaluate_x2*i+1/f*(evaluate_x2*i*Math.sin(X05)+evaluate_x2*i

20、+1*Math.cos(X05);/*a26*/A2*i+15=-evaluate_x2*i;/進行常數(shù)項的初始化for(i=0;i<4;i+)L2*i=xi0-evaluate_x2*i;L2*i+1=xi1-evaluate_x2*i+1;/A的轉(zhuǎn)置矩陣for(i=0;i<8;i+)for(j=0;j<6;j+)A_Tji=Aij;/實現(xiàn)A與AT相乘int k=0;for(i=0;i<6;i+)for(j=0;j<6;j+)AT_Aij=0;for(i=0;i<6;i+)for(k=0;k<6;k+)for(j=0;j<8;j+)AT_Aik

21、+=A_Tij*Ajk;/得到AT*A的逆矩陣存放在inverseAT_A中,此處采用面向?qū)ο蟮姆椒╥nverseAT_A=hengFind.inverse(AT_A,6);/實現(xiàn)矩陣AT_A66與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+=inverseAT_Aij*A_Tjk;/實現(xiàn)result168與l8的相乘，得到結(jié)果放在result26中for(i=0;i&l

22、t;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;System.out.println("n"+"進行第"+num+"次迭代得到Xs,Ys,Zs, ，改正數(shù)分別為：");for(i=0;i<6;i+)System.out.print(result2i+"t");while(Math.abs(result23*206265.0)&

23、gt;6|Math.abs(result24*206265.0)>6|Math.abs(result25*206265.0)>6);System.out.print("n"+"滿足限差條件結(jié)束循環(huán)，最終結(jié)果為："+"n");System.out.println("Xs"+"t"+"Ys"+"t"+"Zs"+"t"+""+"t"+""+"

24、t"+"k");for(i=0;i<6;i+)System.out.print(X0i+"t");/之前代碼所求為外方位元素的六個值/以下代碼所求為單位權(quán)中誤差for(i=0;i<8;i+)result_vi=0;for(i=0;i<8;i+)/這兒可能會出錯for(j=0;j<6;j+)result_vi+=Aij*result2j;for(int q=0;q<6;q+)result_vq-=Lq;for(int m=0;m<8;m+)resultF+=(result_vm*result_vm);resul

25、tF=Math.sqrt(resultF/2);for(int p=0;p<6;p+)standardp=(Math.sqrt(inverseAT_App)*resultF;System.out.println("n"+"各所求量的中誤差為");System.out.print("Xs"+"t"+"Ys"+"t"+"Zs"+"t"+""+"t"+""+"t&qu

26、ot;+"k"+"n");System.out.print(standard0+"t"+standard1+"t"+standard2+"t"+standard3+"t"+standard4+"t"+standard5+"n");以下第二部分為計算逆矩陣的面向?qū)ο蟪绦騪ublic class hengpublic double inverse(double a,int n)double b=new doublen2*n;double c=new doublenn;double det1,yinzhi,bb;int i,j,k,u;for(i=0;i<n;i+)for(j=0;j<2*n;j+)bij=0;for(i=0;i<n;i+)for(j=0;j<n;j+)bij=aij;for(j=0;j