近景摄影测量光束法平差报告(共14页).doc

《近景摄影测量光束法平差报告(共14页).doc》由会员分享，可在线阅读，更多相关《近景摄影测量光束法平差报告(共14页).doc（14页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上近景摄影测量光束法平差报告2011 年 6 月 4 日1 作业目的 -32 外业控制点的观测与解算 -33 近景影像获取 -44 LPS刺点点位 -45 光束法平差与精度评定 -56 总结 -111 作业目的以近景摄影测量大实习为基础，对所摄取近景相片解析处理，以外业控制点的解算成果以及内业LPS平差结果为依据，编写光束法平差程序，由22个控制点的像素坐标及5个“已知控制点”的三维坐标求解其余17个控制点的三维坐标，并评定精度。2 作业条件及数据点号像素坐标列（J）像素坐标行（I）XYZ左片：2650.9892114.93497.4532353.7473299.895

2、382792.4912259.531508.8008342.3524298.6832102791.483740.514508.8138342.3548307.0717163928.5592120.49520.2969353.7531300.1146214890.5842130.45527.9857353.5821300.10371648.6242765.5820003660.4521441.4110004728.563816.58500051965.8952557.99600061910.1051210.0700072767.4553044.53100092774.0591493.061000

3、123319.0112665.417000133312.2861986.582000143298.4681284.901000154055.0522705.029000173808.9851539.018000183715.006962.032000193836.444706.426000204883.392691.6510002247541681000234825.4091018.545000右片：2670.9482129.967497.4532353.7473299.895382346.4432264.542508.8008342.3524298.6832102361.448691.079

4、508.8138342.3548307.0717164088.4192115.427520.2969353.7531300.1146205203.4412736.112527.9857353.5821300.10371666.1032764.8820003685.4031472.5740004754.414860.65600051652.4312568.50300061600.0581207.01400072312.0273077.96400092334.4721470.473000123083.1932691.257000133083.0661976.987000143072.4281240

5、.419000154230.5272741.493000173956.4451498.102000183852.033888.02000194075.943613.315000215214.4922119.571000225144.4631507.538000235139.982903.989000外方位元素初始值：(左)Xs1 = 497.9149,Ys1 = 301.2754,Zs1 = 297.2430,Q1=14.9560,W1=4.7765,K1=-0.0308,(右)Xs2 = 509.6501,Ys2 = 301.3448,Zs2 = 297.4727,Q2=2.1777 ,W2

6、=4.5969,K2=0.0791,相片主距：50mm3 平差思想31 共线条件方程本次作业中，光束法平差基于如下共线条件方程：该共线方程是描述摄影中心S、像点a以及物点A位于一直线上的关系式。像点a在成像过程中存在某种系统误差，其改正数（x，y）添加在上式左边，并有：将共线条件方程线性化，其中：32 像点坐标误差方程式由于本作业采用光束法平差，两张照片共44个点（每张22个），每个点可列2个像点坐标误差方程（vx，vy），故总共88个误差方程。而每个误差方程中：有12个外方位元素改正（每张相片6个），6个内方位元素改正（每张相片3个），51个加密点坐标改正（5个控制点坐标改正为0，故只有有1

7、7个加密点的三维坐标改正：dX，dY,dZ），4个畸变参数改正（每张相片2个，k1，k2）。33 平差由88个像点坐标误差方程式组成方程组：,其中：V阵为像点坐标误差改正，88行*1列；A阵为系数阵，88行*73列；X阵伟未知数阵，73行*1列（未知数包括12个外方位元素、6个内方位元素、51个加密点三维坐标、4个畸变参数）L阵为常数项阵，88行*1列由式X=（AtA）-1 ATL解求未知数即可4 程序#include #include #include #include#includeconst int N=90;/求转置矩阵templatevoid Transpose(T1*mat1,T2

8、*mat2,int a,int b)int i,j;for(i=0;ia;i+)for(j=0;jb;j+)mat2ji=mat1ij;return;/求矩阵的乘积templatevoid Array_mul(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resultij=0;for(k=0;kb;k+)resultij+=mat1ik*mat2kj;return;/求逆矩阵inverse(double ANN,int m)int i=0,j=0,k=0;double

9、 C2020,B2020;for(i=0;i2*m;i+) for(j=0;j2*m;j+) if(i=j) Cij=1.0; else Cij=0.0;for(i=0;im;i+) for(j=0;jm;j+)Bij=Aij;for(i=0;im;i+) for(j=m;j2*m;j+)Bij=Cij-m;cout.precision(5);for(k=0;km;k+)for(i=k;i=0;j-)if(Bik!=0)Bij=Bij/Bik; for(i=m-1;ik;i-)if(Bik=0)continue; for(j=0;j2*m;j+) Bij=Bij-Bkj; for(k=1;km

10、;k+) for(i=0;i=i;j-) Bij=Bij-Bik*Bkj;for(i=0;im;i+) for(j=0;jm;j+)Aji=Bjm+i; return 1; void main()/定义两张相片共个控制点和加密点的像素坐标（J1,I1,J2,I2）和像平面直角坐标（x1,y1,x2,y2），及地面坐标（X,Y,Z）double J1N=0.0,I1N=0.0,J2N=0.0,I2N=0.0;double x1N=0.0,x2N=0.0,z1N=0.0,z2N=0.0;double XN=0.0,YN=0.0,ZN=0.0; /导入控制点坐标数据ifstream infile;i

11、nfile.open(控制点坐标数据.txt); if(infile.is_open() while(!infile.eof ()for(int i=0;iJ1i;infile.ignore(1); infileI1i;infile.ignore(1); infileXi;infile.ignore(1); infileYi;infile.ignore(1); infileZi;infile.ignore(1); infile.close(); coutJ11 J143endl;/像素坐标转化为直角坐标for (int j = 0; j 44; j+)x1j = (J1j - 5616/2)

12、* 6. / 1000 ;z1j = (3744/2 - I1j) * 6. / 1000 ;coutx11 x143endl;/左右相片外方位元素初始值及摄影机主距double Xs1,Ys1,Zs1,Q1,W1,K1,f1,x01,z01,k11,k21;double Xs2,Ys2,Zs2,Q2,W2,K2,f2,x02,z02,k12,k22;Xs1 = 497.9149,Ys1 = 301.2754,Zs1 = 297.2430,Q1=14.9560,W1=4.7765,K1=-0.0308,f1 = 50, x01 = 0, z01 = 0, k11 = 0, k21 = 0;Xs

13、2 = 509.6501,Ys2 = 301.3448,Zs2 = 297.4727,Q2=2.1777 ,W2=4.5969,K2=0.0791,f2 = 50, x02 = 0, z02 = 0, k12 = 0, k22 = 0;double rrN=0.0,dxN=0.0,dzN=0.0;double XX1,YY1,ZZ1,XX2,YY2,ZZ2;/定义旋转矩阵，系数矩阵,常数项和改正数double R133=0.0,R233=0.0,A1NN=0.0,l1NN=0.0,dNN=0.0;int t=0;coutXs1endl;/组成旋转矩阵 R100=cos(Q1)*cos(K1)-

14、sin(Q1)*sin(W1)*sin(K1); R101=cos(W1)*sin(Q1); R102=-cos(Q1)*sin(K1)-sin(Q1)*sin(W1)*cos(K1); R110=-sin(Q1)*cos(K1)-cos(Q1)*sin(W1)*sin(K1); R111=cos(Q1)*cos(W1); R112=sin(Q1)*sin(K1)-cos(Q1)*sin(W1)*cos(K1); R120=cos(W1)*sin(K1); R121=-sin(W1); R122=cos(W1)*cos(K1); R200=cos(Q2)*cos(K2)-sin(Q2)*sin

15、(W2)*sin(K2); R201=cos(W2)*sin(Q2); R202=-cos(Q2)*sin(K2)-sin(Q2)*sin(W2)*cos(K2); R210=-sin(Q2)*cos(K2)-cos(Q2)*sin(W2)*sin(K2); R211=cos(Q2)*cos(W2); R212=sin(Q2)*sin(K2)-cos(Q2)*sin(W2)*cos(K2); R220=cos(W2)*sin(K2); R221=-sin(W2); R222=cos(W2)*cos(K2); for(int k=0;k4;k+) XX1=R100*(Xk-Xs1)+R110*(

16、Yk-Ys1)+R120*(Zk-Zs1); YY1=R101*(Xk-Xs1)+R111*(Yk-Ys1)+R121*(Zk-Zs1); ZZ1=R102*(Xk-Xs1)+R112*(Yk-Ys1)+R122*(Zk-Zs1); XX2=R200*(Xk-Xs1)+R210*(Yk-Ys1)+R220*(Zk-Zs1); YY2=R201*(Xk-Xs1)+R211*(Yk-Ys1)+R221*(Zk-Zs1); ZZ2=R202*(Xk-Xs1)+R212*(Yk-Ys1)+R222*(Zk-Zs1); for(int p=0;p22;p+) rrp=(x1p-x01)*(x1p-x01

17、)+(z1p-z01)*(z1p-z01); dxp=(x1p-x01)*(k11*(x1p-x01)*(x1p-x01)+k21*(x1p-x01)*(x1p-x01)*(x1p-x01)*(x1p-x01); dzp=(z1p-z01)*(k11*(z1p-z01)*(z1p-z01)+k21*(z1p-z01)*(z1p-z01)*(z1p-z01)*(z1p-z01); coutrr1endl; for(int q=22;q44;q+) rrq=(x1q-x02)*(x1q-x02)+(z1q-z02)*(z1q-z02); dxq=(x1q-x02)*(k12*(x1q-x02)*(

18、x1q-x02)+k22*(x1q-x02)*(x1q-x02)*(x1q-x02)*(x1q-x02); dzq=(z1q-z02)*(k12*(z1q-z02)*(z1q-z02)+k22*(z1q-z02)*(z1q-z02)*(z1q-z02)*(z1q-z02); coutrr22endl; /计算系数阵(88*73)和常数项 /左片 for(int i=0,H=0;i22;i+) l1H0=x1i-x01+f1*XX1/YY1-dxi;l1H+10=z1i-z01+f1*ZZ1/YY1-dzi;A1H0=(R101*(x1i-x01)-R100*f1)/YY1;/x/Xs=-x/X

19、A1H1=(R111*(x1i-x01)-R110*f1)/YY1;/x/YsA1H2=(R121*(x1i-x01)-R120*f1)/YY1;/x/ZsA1H3=(z1i-z01)*sin(W1)-(x1i-x01)*(x1i-x01)*cos(K1)-(z1i-z01)*sin(K1)/f1+f1*cos(K1)*cos(W1);/x/QA1H4=-f1*sin(K1)-(x1i-x01)*(x1i-x01)*sin(K1)+(z1i-z01)*cos(K1)/f1;/x/WA1H5=z1i-z01;/x/KA1H6=(x1i-x01)/f1;/x/fA1H7=1;/x/x0A1H8=0

20、;/x/z0A1H9=(x1i-x01)*rri;/x/k1A1H10=(x1i-x01)*rri*rri;/x/k2A1H11=(R201*(x1i-x02)-R200*f2)/YY2;/x/Xs=-x/XA1H12=(R211*(x1i-x02)-R210*f2)/YY2;/x/YsA1H13=(R221*(x1i-x02)-R220*f2)/YY2;/x/ZsA1H14=(z1i-z02)*sin(W2)-(x1i-x02)*(x1i-x02)*cos(K2)-(z1i-z02)*sin(K2)/f2+f2*cos(K2)*cos(W2);/x/QA1H15=-f2*sin(K2)-(x

21、1i-x02)*(x1i-x02)*sin(K2)+(z1i-z02)*cos(K2)/f2;/x/WA1H16=z1i-z02;/x/KA1H17=(x1i-x02)/f2;/x/fA1H18=1;/x/x0A1H19=0;/x/z0A1H20=(x1i-x02)*rri;/x/k1A1H21=(x1i-x02)*rri*rri;/x/k2A1H+10=(R101*(z1i-z01)-R102*f1)/YY1;/z/Xs=-z/XA1H+11=(R111*(z1i-z01)-R112*f1)/YY1;/z/Ys=-z/YA1H+12=(R121*(z1i-z01)-R122*f1)/YY1;

22、/z/Zs=-z/ZA1H+13=-(x1i-x01)*sin(W1)-(z1i-z01)*(x1i-x01)*cos(K1)-(z1i-z01)*sin(K1)/f1-f1*sin(K1)*cos(W1);/z/QA1H+14=-f1*cos(K1)-(z1i-z01)*(x1i-x01)*sin(K1)+(z1i-z01)*cos(K1)/f1;/z/WA1H+15=-(x1i-x01);/z/KA1H+16=(z1i-z01)/f1;/z/fA1H+17=0;/z/x0A1H+18=1;/z/z0A1H+19=(z1i-z01)*rri;/x/k1A1H+110=(z1i-z01)*rr

23、i*rri;/x/k2A1H+111=(R201*(z1i-z02)-R202*f2)/YY2;/z/Xs=-z/XA1H+112=(R211*(z1i-z02)-R212*f2)/YY2;/z/Ys=-z/YA1H+113=(R221*(z1i-z02)-R222*f2)/YY2;/z/Zs=-z/ZA1H+114=-(x1i-x02)*sin(W2)-(z1i-z02)*(x1i-x02)*cos(K2)-(z1i-z02)*sin(K2)/f2-f2*sin(K2)*cos(W2);/z/QA1H+115=-f2*cos(K2)-(z1i-z02)*(x1i-x02)*sin(K2)+(

24、z1i-z02)*cos(K2)/f2;/z/WA1H+116=-(x1i-x02);/z/KA1H+117=(z1i-z02)/f2;/z/fA1H+118=0;/z/x0A1H+119=1;/z/z0A1H+120=(z1i-z02)*rri;/x/k1A1H+121=(z1i-z02)*rri*rri;/x/k2H=H+2; /右片 for(int ii=22,HH=44;ii44;ii+) l1HH0=x1ii-x01+f1*XX1/YY1-dxii;l1HH+10=z1ii-z01+f1*ZZ1/YY1-dzii;A1HH0=(R101*(x1ii-x01)-R100*f1)/YY1

25、;/x/Xs=-x/XA1HH1=(R111*(x1ii-x01)-R110*f1)/YY1;/x/YsA1HH2=(R121*(x1ii-x01)-R120*f1)/YY1;/x/ZsA1HH3=(z1ii-z01)*sin(W1)-(x1ii-x01)*(x1ii-x01)*cos(K1)-(z1ii-z01)*sin(K1)/f1+f1*cos(K1)*cos(W1);/x/QA1HH4=-f1*sin(K1)-(x1ii-x01)*(x1ii-x01)*sin(K1)+(z1ii-z01)*cos(K1)/f1;/x/WA1HH5=z1ii-z01;/x/KA1HH6=(x1ii-x0

26、1)/f1;/x/fA1HH7=1;/x/x0A1HH8=0;/x/z0A1HH9=(x1ii-x01)*rrii;/x/k1A1HH10=(x1ii-x01)*rrii*rrii;/x/k2A1HH11=(R201*(x1ii-x02)-R200*f2)/YY2;/x/Xs=-x/XA1HH12=(R211*(x1ii-x02)-R210*f2)/YY2;/x/YsA1HH13=(R221*(x1ii-x02)-R220*f2)/YY2;/x/ZsA1HH14=(z1ii-z02)*sin(W2)-(x1ii-x02)*(x1ii-x02)*cos(K2)-(z1ii-z02)*sin(K2

27、)/f2+f2*cos(K2)*cos(W2);/x/QA1HH15=-f2*sin(K2)-(x1ii-x02)*(x1ii-x02)*sin(K2)+(z1ii-z02)*cos(K2)/f2;/x/WA1HH16=z1ii-z02;/x/KA1HH17=(x1ii-x02)/f2;/x/fA1HH18=1;/x/x0A1HH19=0;/x/z0A1HH20=(x1ii-x02)*rrii;/x/k1A1HH21=(x1ii-x02)*rrii*rrii;/x/k2A1HH+10=(R101*(z1ii-z01)-R102*f1)/YY1;/z/Xs=-z/XA1HH+11=(R111*(

28、z1ii-z01)-R112*f1)/YY1;/z/Ys=-z/YA1HH+12=(R121*(z1ii-z01)-R122*f1)/YY1;/z/Zs=-z/ZA1HH+13=-(x1ii-x01)*sin(W1)-(z1ii-z01)*(x1ii-x01)*cos(K1)-(z1ii-z01)*sin(K1)/f1-f1*sin(K1)*cos(W1);/z/QA1HH+14=-f1*cos(K1)-(z1ii-z01)*(x1ii-x01)*sin(K1)+(z1ii-z01)*cos(K1)/f1;/z/WA1HH+15=-(x1ii-x01);/z/KA1HH+16=(z1ii-z0

29、1)/f1;/z/fA1HH+17=0;/z/x0A1HH+18=1;/z/z0A1HH+19=(z1ii-z01)*rrii;/x/k1A1HH+110=(z1ii-z01)*rrii*rrii;/x/k2A1HH+111=(R201*(z1ii-z02)-R202*f2)/YY2;/z/Xs=-z/XA1HH+112=(R211*(z1ii-z02)-R212*f2)/YY2;/z/Ys=-z/YA1HH+113=(R221*(z1ii-z02)-R222*f2)/YY2;/z/Zs=-z/ZA1HH+114=-(x1ii-x02)*sin(W2)-(z1ii-z02)*(x1ii-x02

30、)*cos(K2)-(z1ii-z02)*sin(K2)/f2-f2*sin(K2)*cos(W2);/z/QA1HH+115=-f2*cos(K2)-(z1ii-z02)*(x1ii-x02)*sin(K2)+(z1ii-z02)*cos(K2)/f2;/z/WA1HH+116=-(x1ii-x02);/z/KA1HH+117=(z1ii-z02)/f2;/z/fA1HH+118=0;/z/x0A1HH+119=1;/z/z0A1HH+120=(z1ii-z02)*rrii;/x/k1A1HH+121=(z1ii-z02)*rrii*rrii;/x/k2HH=HH+2; /计算左片外方位元素

31、Xs1,Ys1,Zs1,Q,W,Kdouble A1TNN=0,A1TA1NN=0,A1TL1NN=0;Transpose(A1,A1T,N,N);Array_mul(A1T,A1,A1TA1,N,N,N); inverse(A1TA1,N);Array_mul(A1T,l1,A1TL1,N,N,N);Array_mul(A1TA1,A1TL1,d,N,N,N);Xs1=Xs1+d00; Ys1=Ys1+d10; Zs1=Zs1+d20;Q1=Q1+d30; W1=W1+d40; K1=K1+d50;f1=f1+d60; x01=x01+d70; z01=z01+d80; k11=k11+d9

32、0; k21=k21+d100;Xs2=Xs2+d110; Ys2=Ys1+d120; Zs2=Zs2+d130;Q2=Q2+d140; W2=W2+d150; K2=K2+d160;f2=f2+d170; x02=x02+d180; z02=z02+d190;k12=k11+d90; k22=k22+d100;cout迭代次数：tendl;cout外方位元素的直线元素为endl;cout Xs1 = Xs1 , Xs2 = Xs2endl;cout Ys1 = Ys1 , Ys2 = Ys2endl;cout Zs1 = Zs1 , Zs2 = Zs2endl;cout Q1 = Q1 ,

33、Q2 = Q2endl;cout W1 = Q1 , W2 = W2endl; cout K1 = K1 , K2 = K2endl;coutendl;cout f1 = f1 , f2 = f2endl;cout x01 = f1 , x02 = f2endl;cout z01 = f1 , z02 = f2endl;coutendl;cout k11 = k11 , k12 = k12endl;cout k21 = k21 , k22 = k22endl;5 结果在这里要请求老师原谅，很抱歉在编制程序过程中，出现了陷入死循环的问题，一直修改调试也未能解决问题，导致无法计算出结果。整个程序也

34、还很不完整，当然，我知道自己的程序还存在很多问题，但是由于时间问题，未能及时一一解决，还望老师原谅！6 总结本次程序编制作业完成得很失败，由于一直没有很多精力和时间来调试和修改程序，最后无法得出结果，再次表示抱歉。但是整个平差的过程我已经掌握，当然从懵懂到完全了解掌握也经历了不少波折，开始的时候还没有弄清楚“光束法”的含义，是应该两张照片同时解算，而一开始我混淆概念，用了后方交会的方法求解，后来才发现这样只能解算单张相片的外方位元素。由于本次近景摄影测量测量的为6教墙面上的点，相片平面位于XZ面而非XY面，故共线条件方程也有所改变，这点也是有同学发现后才察觉到。根据老师修改后的共线条件方程及相关参数，开始列像点坐标误差方程，一开始没有想象到误差方程会这么庞大，但也是在和同学的讨论中最后得出了正确的误差方程及系数阵等的列法。本次作业原理并不难理解，难点在“光束法”及列误差方程和系数阵上吧，平差过程也不复杂，不过数据量较庞大。通过这次作业，我很深刻地理解了光束法平差的思想和过程，也掌握了它的方法，但是很遗憾没能得到结果，熬过考试周到暑假我一定会抽出时间把它解决掉，完善自己的程序，从中吸取经验教训，不断提升自我。专心-专注-专业