近景摄影测量光束法平差报告

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.895382792.491225

2、9.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.061000123319.011266

3、5.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.079508.8138342.3

4、548307.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.419000154230

5、.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=4.5969,K2=0.

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

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 <iomanip.h>#include <stdlib.h>#include <math.h>#include<fstream.h>#include<

8、;iostream.h>const int N=90;/求转置矩阵template<typename T1,typename T2>void Transpose(T1*mat1,T2*mat2,int a,int b)int i,j;for(i=0;i<a;i+)for(j=0;j<b;j+)mat2ji=mat1ij;return;/求矩阵的乘积template<typename T1,typename T2>void Array_mul(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j

9、,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;/求逆矩阵inverse(double ANN,int m)int i=0,j=0,k=0;double C2020,B2020;for(i=0;i<2*m;i+) for(j=0;j<2*m;j+) if(i=j) Cij=1.0; else Cij=0.0;for(i=0;i<m;i+) for(j=0;j<m;j+)Bij=Aij;for(i=0;i<m;i+)

10、for(j=m;j<2*m;j+)Bij=Cij-m;cout.precision(5);for(k=0;k<m;k+)for(i=k;i<m;i+)for(j=2*m-1;j>=0;j-)if(Bik!=0)Bij=Bij/Bik; for(i=m-1;i>k;i-)if(Bik=0)continue; for(j=0;j<2*m;j+) Bij=Bij-Bkj; for(k=1;k<m;k+) for(i=0;i<k;i+) for(j=2*m-1;j>=i;j-) Bij=Bij-Bik*Bkj;for(i=0;i<m;i+)

11、for(j=0;j<m;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;infile.open("控制点坐标数据.txt"); if(infile.i

12、s_open() while(!infile.eof ()for(int i=0;i<44;i+) infile>>J1i;infile.ignore(1); infile>>I1i;infile.ignore(1); infile>>Xi;infile.ignore(1); infile>>Yi;infile.ignore(1); infile>>Zi;infile.ignore(1); infile.close(); cout<<J11<<" "<<J143<&l

13、t;endl;/像素坐标转化为直角坐标for (int j = 0; j < 44; j+)x1j = (J1j - 5616/2) * 6.410256 / 1000 ;z1j = (3744/2 - I1j) * 6.410256 / 1000 ;cout<<x11<<" "<<x143<<endl;/左右相片外方位元素初始值及摄影机主距double Xs1,Ys1,Zs1,Q1,W1,K1,f1,x01,z01,k11,k21;double Xs2,Ys2,Zs2,Q2,W2,K2,f2,x02,z02,k12,k

14、22;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;Xs2 = 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,YY

15、2,ZZ2;/定义旋转矩阵，系数矩阵,常数项和改正数double R133=0.0,R233=0.0,A1NN=0.0,l1NN=0.0,dNN=0.0;int t=0;cout<<Xs1<<endl;/组成旋转矩阵 R100=cos(Q1)*cos(K1)-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(

16、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(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=s

17、in(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;k<4;k+) XX1=R100*(Xk-Xs1)+R110*(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*

18、(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;p<22;p+) rrp=(x1p-x01)*(x1p-x01)+(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*

19、(z1p-z01)*(z1p-z01)*(z1p-z01)*(z1p-z01); cout<<rr1<<endl; for(int q=22;q<44;q+) rrq=(x1q-x02)*(x1q-x02)+(z1q-z02)*(z1q-z02); dxq=(x1q-x02)*(k12*(x1q-x02)*(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)*(z1

20、q-z02); cout<<rr22<<endl; /计算系数阵(88*73)和常数项 /左片 for(int i=0,H=0;i<22;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/XA1H1=(R111*(x1i-x01)-R110*f1)/YY1;/x/YsA1H2=(R121*(x1i-x01)-R120*f1)/YY1;/x/ZsA1H3=(z1i-z01)*sin(W1)-(x1i-x01)

21、*(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;/x/z0A1H9=(x1i-x01)*rri;/x/k1A1H10=(x1i-x01)*rri*rri;/x/k2A1H11=(R201*(x1i-x02)-R200*f2)/YY2;/x/Xs=-x/XA1H12=

22、(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)-(x1i-x02)*(x1i-x02)*sin(K2)+(z1i-z02)*cos(K2)/f2;/x/WA1H16=z1i-z02;/x/KA1H17=(x1i-x02)/f2;/x/fA1H18=1;/x/x0A1H19

23、=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;/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/QA1

24、H+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)*rri*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+1

25、13=(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)+(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=

26、(z1i-z02)*rri;/x/k1A1H+121=(z1i-z02)*rri*rri;/x/k2H=H+2; /右片 for(int ii=22,HH=44;ii<44;ii+) l1HH0=x1ii-x01+f1*XX1/YY1-dxii;l1HH+10=z1ii-z01+f1*ZZ1/YY1-dzii;A1HH0=(R101*(x1ii-x01)-R100*f1)/YY1;/x/Xs=-x/XA1HH1=(R111*(x1ii-x01)-R110*f1)/YY1;/x/YsA1HH2=(R121*(x1ii-x01)-R120*f1)/YY1;/x/ZsA1HH3=(z1ii-z

27、01)*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-x01)/f1;/x/fA1HH7=1;/x/x0A1HH8=0;/x/z0A1HH9=(x1ii-x01)*rrii;/x/k1A1HH10=(x1ii-x01)*rrii*rrii;/x/k2A1HH11=(R2

28、01*(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)/f2+f2*cos(K2)*cos(W2);/x/QA1HH15=-f2*sin(K2)-(x1ii-x02)*(x1ii-x02)*sin(K2)+(z1ii-z02)*cos(K2)/f2;/x/WA1H

29、H16=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*(z1ii-z01)-R112*f1)/YY1;/z/Ys=-z/YA1HH+12=(R121*(z1ii-z01)-R122*f1)/YY1;/z/Zs=-z/ZA1HH+13=-(x1ii-x01)*sin(W

30、1)-(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-z01)/f1;/z/fA1HH+17=0;/z/x0A1HH+18=1;/z/z0A1HH+19=(z1ii-z01)*rrii;/x/k1A1HH+110=(z1ii-z01)*rrii*rrii;/x/k2A1

31、HH+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)*cos(K2)-(z1ii-z02)*sin(K2)/f2-f2*sin(K2)*cos(W2);/z/QA1HH+115=-f2*cos(K2)-(z1ii-z02)*(x1ii-x02)*sin(K2)+

32、(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; /计算左片外方位元素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); i

33、nverse(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+d90; k21=k21+d100;Xs2=Xs2+d110; Ys2=Ys1+d120; Zs2=Zs2+d130;Q2=Q2+d140; W2=W2+d150; K2=K2+d160;f2=f2+d170; x0

34、2=x02+d180; z02=z02+d190;k12=k11+d90; k22=k22+d100;cout<<"迭代次数："<<t<<endl;cout<<"外方位元素的直线元素为"<<endl;cout<<" Xs1 = "<<Xs1<<" , Xs2 = "<<Xs2<<endl;cout<<" Ys1 = "<<Ys1<<"

35、; , Ys2 = "<<Ys2<<endl;cout<<" Zs1 = "<<Zs1<<" , Zs2 = "<<Zs2<<endl;cout<<" Q1 = "<<Q1<<" , Q2 = "<<Q2<<endl;cout<<" W1 = "<<Q1<<" , W2 = "<<W2<<endl; cou