单像空间后方交会程序《Java版》

1、本科实验报告 课程名称：单像空间后方交会编程实现 姓 名：李恒松 学 院：地学学院 专 业：测绘工程一班 学 号： 指导老师：陈强 2015年 5 月 20 日目录一、作业任务- 3 -二、计算原理- 3 -三、算法流程- 7 -四、源程序- 8 -五、计算结果- 8 -六、结果分析- 8 -七、心得与体会- 8 -八、附页- 8 -Java程序- 8 -一、 作业任务已知条件：摄影机主距f=153.24mm，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。 点号像点坐标地面坐标x(mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.993

2、6589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31以单像空间后方交会方法，求解该像片的外方位元素。二、 计算原理1. 获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄测坐标。2. 测量控制点的像点坐标并作系统误差改正。3. 确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即式中：为摄影比例尺分母；为控制点个数 4. 用

3、三个角元素的初始值按下式计算各方向余弦值，组成旋转矩阵其中：5. 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标，带入以下共线方程式， 逐点计算像点坐标的近似值、，将题目所给数据代入上式可得像点坐标的近似数据。6. 逐点计算误差方程式的系数和常数项，组成误差方程式。其中其他项类似可得。由常数项计算公式：将题目所给数据代入上式，可得常数项。7. 计算法方程的是系数矩阵和常数项，组成法方程式。8. 解法方程，求得外方位元素的改正数。9. 用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。式中：代表迭代次数。10. 将求得的外方位元素改正数与规定的限差比较，若小于限差

4、，则迭代结束。否则用新的近似值重复49满足要求为止。11. 由误差方程的计算式：以及单位权中误差的计算式（式中：表示多余观测数）以及平差中6个参数的协因数阵最终得到参数的中误差为一、 算法流程否是否是输入原始数据像点坐标计算，系统误差改正确定外方位元素初始值组成旋转矩阵R逐点组成误差方程式并法化所有像点完否计算改正后的外方位元素计算中误差，输出成果，结束解法方程，求外方位元素改正数外方位元素改正数是否小于限差是结束并提示错误信息迭代次数小于n二、 源程序源程序代码请见附页三、 计算结果在经过三次迭代之后得到的最终成果如表1所示表1，最终计算结果Xs（米） Ys（米）Zs（米）（弧度）（弧度）（

5、弧度）39795.27476.47572.6-0.0.-0.四、 结果分析由计算结果可知在拍摄照片瞬间，摄影中心在地面摄影测量坐标系中的坐标为（39795.，27476.，7572.）（单位：米），航向倾角为-0.弧度，旁向倾角为0.弧度，像片旋角为0.弧度。表2，精度评定结果参数中误差Xs1.957Ys1.7318Zs0.769861.19645E-41.39593E-47.2504E-5五、 心得与体会通过这次大作业，我的心得体会以下几点：第一：编程确实是一件非诚艰难的事情，即使整个单像空间后方交会实验的原理非常简单，清晰明了，但是要用编程把它表示出来还是一件非常艰难的事情。我遇到的最大的

6、难点是在矩阵求逆的步骤上，为了解决这个问题。我用了两天的时间，还请教了数学学院的同学，程序里求逆矩阵的方法是来自那位同学的。这次编程也让我对于线性代数中的求逆矩阵的问题理解的更加深刻。第二：在发现问题，解决问题的过程中，我也认识到了自己基础知识的不足，对于Java语言的不熟悉让我在很多常见的bug面前无能为力。以后要加强编程语言的学习。第三：很多东西看似简单，其实细节处才体现出能力。把小的细节处理好，才能有大的作品。六、 附页1. Java程序以下第一部分为主程序import java.util.*;import java.math.*;public class li public stati

7、c void main(String args)heng hengFind=new heng();double x=-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.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=

8、1/40000.0;double R=new double33;double evaluate_x=new double8;double A=new double86;double A_T=new double68;double L=new double8;double pi=3.97932;double AT_A=new double66;double inverseAT_A=new double66;/inverseAT_A=hengFind.inverse(AT_A,6);double result_v=new double8;/用来存放用于计算单位权中误差的v值double stand

9、ard=new double6;double resultF=0;/用来存放v的转置矩阵double result1=new double68;double result2=new double6;double sumXYZ=new double3;System.out.println(已知像点坐标为： );for(i=0;i4;i+)for(j=0;j2;j+)if(j=0)System.out.println(x+(i+1)+= +t+xij);elseSystem.out.println(y+(i+1)+= +t+xij);System.out.println(已知地面四个点的坐标为：)

10、;for(i=0;i4;i+)for(j=0;j3;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;j3;j+)for(i=0;i4;i+)sumXYZj+=Xij;for(i=0;i2;i+)X0i=sumXYZi/4;X0i=1/a*f+sumXYZ2/4.0;doR00=Math.cos(X03)*Math.cos(X05)-Math.sin(X03)

11、*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.sin(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=-M

12、ath.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)/(R02*(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-X0

13、2);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)/(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

14、-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-X02)/(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*(

15、X32-X02);for(i=0;i4;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)+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*evalu

16、ate_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=(R21*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*

17、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)+evaluate_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(X

18、05)-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+1*Math.cos(X05);/*a26*/A2*i+15=-evaluate_x2*i;/进行常数项的初始化for(i=0;i4;i+)L2*i=xi0-evaluate_x2*i;L2*i+1=xi1-evaluate_x2*i+1;/A的转置矩阵for(i=0;i8;i+)

19、for(j=0;j6;j+)A_Tji=Aij;/实现A与AT相乘int k=0;for(i=0;i6;i+)for(j=0;j6;j+)AT_Aij=0;for(i=0;i6;i+)for(k=0;k6;k+)for(j=0;j8;j+)AT_Aik+=A_Tij*Ajk;/得到AT*A的逆矩阵存放在inverseAT_A中,此处采用面向对象的方法inverseAT_A=hengFind.inverse(AT_A,6);/实现矩阵AT_A66与AT68的相乘,结果存放在result168中for(i=0;i6;i+)for(j=0;j8;j+)result1ij=0;for(i=0;i6;i

20、+)for(k=0;k8;k+)for(j=0;j6;j+)result1ik+=inverseAT_Aij*A_Tjk;/实现result168与l8的相乘，得到结果放在result26中for(i=0;i6;i+)result2i=0;for(i=0;i6;i+)for(j=0;j8;j+)result2i+=result1ij*Lj;num+;for(i=0;i6;i+)X0i=X0i+result2i;System.out.println(n+进行第+num+次迭代得到Xs,Ys,Zs, ，改正数分别为：);for(i=0;i6|Math.abs(result24*.0)6|Math.

21、abs(result25*.0)6);System.out.print(n+满足限差条件结束循环，最终结果为：+n);System.out.println(Xs+t+Ys+t+Zs+t+t+t+k);for(i=0;i6;i+)System.out.print(X0i+t);/之前代码所求为外方位元素的六个值/以下代码所求为单位权中误差for(i=0;i8;i+)result_vi=0;for(i=0;i8;i+)/这儿可能会出错for(j=0;j6;j+)result_vi+=Aij*result2j;for(int q=0;q6;q+)result_vq-=Lq;for(int m=0;m

22、8;m+)resultF+=(result_vm*result_vm);resultF=Math.sqrt(resultF/2);for(int p=0;p6;p+)standardp=(Math.sqrt(inverseAT_App)*resultF;System.out.println(n+各所求量的中误差为);System.out.print(Xs+t+Ys+t+Zs+t+t+t+k+n);System.out.print(standard0+t+standard1+t+standard2+t+standard3+t+standard4+t+standard5+n);以下第二部分为计算逆矩阵的面向对象程序public 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;in;i+)for(j=0;j2*n;j+)bij=0;for(i=0;in;i+)for(j=0;jn;j+)bij=aij;for(j=0;jn;j+)