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

近景摄影测量光束法平差报告2011年6月4日1作业目的------------------------------------------------------------------------------------32外业控制点的观测与解算--------------------------------------------------------33近景影像获取----------------------------------------------------------------------------44LPS刺点点位------------------------------------------------------------------45光束法平差与精度评定------------------------------------------------------------56总结---------------------------------------------------------------------------------------------111作业目的以近景摄影测量大实习为基础，对所摄取近景相片解析处理，以外业控制点的解算成果以及内业LPS平差结果为依据，编写光束法平差程序，由22个控制点的像素坐标及5个“已知控制点”的三维坐标求解其余17个控制点的三维坐标，并评定精度。2作业条件及数据点号像素坐标列（J）像素坐标行（I）XYZ左片：2650.9892114.93497.4532353.7473299.895382792.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.061000123319.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.079508.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.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=4.5969,K2=0.0791,相片主距：50mm3平差思想3．1共线条件方程本次作业中，光束法平差基于如下共线条件方程：该共线方程是描述摄影中心S、像点a以及物点A位于一直线上的关系式。像点a在成像过程中存在某种系统误差，其改正数（△x，△y）添加在上式左边，并有：将共线条件方程线性化，其中：3．2像点坐标误差方程式由于本作业采用光束法平差，两张照片共44个点（每张22个），每个点可列2个像点坐标误差方程（vx，vy），故总共88个误差方程。而每个误差方程中：有12个外方位元素改正（每张相片6个），6个内方位元素改正（每张相片3个），51个加密点坐标改正（5个控制点坐标改正为0，故只有有17个加密点的三维坐标改正：dX，dY,dZ），4个畸变参数改正（每张相片2个，k1，k2）。3．3平差由88个像点坐标误差方程式组成方程组：VAXL,其中：V阵为像点坐标误差改正，88行*1列；A阵为系数阵，88行*73列；X阵伟未知数阵，73行*1列（未知数包括12个外方位元素、6个内方位元素、51个加密点三维坐标、4个畸变参数）L阵为常数项阵，88行*1列由式X=（AtA）-1ATL解求未知数即可4程序#includeiomanip.h#includestdlib.h#includemath.h#includefstream.h#includeiostream.hconstintN=90;//求转置矩阵templatetypenameT1,typenameT2voidTranspose(T1*mat1,T2*mat2,inta,intb){inti,j;for(i=0;ia;i++)for(j=0;jb;j++)mat2[j][i]=mat1[i][j];return;}//求矩阵的乘积templatetypenameT1,typenameT2voidArray_mul(T1*mat1,T2*mat2,T2*result,inta,intb,intc){inti,j,k;for(i=0;ia;i++){for(j=0;jc;j++){result[i][j]=0;for(k=0;kb;k++)result[i][j]+=mat1[i][k]*mat2[k][j];}}return;}//求逆矩阵inverse(doubleA[N][N],intm){inti=0,j=0,k=0;doubleC[20][20],B[20][20];for(i=0;i2*m;i++)for(j=0;j2*m;j++){if(i==j)C[i][j]=1.0;elseC[i][j]=0.0;}for(i=0;im;i++)for(j=0;jm;j++)B[i][j]=A[i][j];for(i=0;im;i++)for(j=m;j2*m;j++)B[i][j]=C[i][j-m];cout.precision(5);for(k=0;km;k++){for(i=k;im;i++)for(j=2*m-1;j=0;j--){if(B[i][k]!=0)B[i][j]=B[i][j]/B[i][k];}for(i=m-1;ik;i--){if(B[i][k]==0)continue;for(j=0;j2*m;j++)B[i][j]=B[i][j]-B[k][j];}}for(k=1;km;k++)for(i=0;ik;i++)for(j=2*m-1;j=i;j--)B[i][j]=B[i][j]-B[i][k]*B[k][j];for(i=0;im;i++)for(j=0;jm;j++)A[j][i]=B[j][m+i];return1;}voidmain(){//定义两张相片共个控制点和加密点的像素坐标（J1,I1,J2,I2）和像平面直角坐标（x1,y1,x2,y2），及地面坐标（X,Y,Z）doubleJ1[N]={0.0},I1[N]={0.0},J2[N]={0.0},I2[N]={0.0};doublex1[N]={0.0},x2[N]={0.0},z1[N]={0.0},z2[N]={0.0};doubleX[N]={0.0},Y[N]={0.0},Z[N]={0.0};//导入控制点坐标数据ifstreaminfile;infile.open(控制点坐标数据.txt);if(infile.is_open()){while(!infile.eof()){for(inti=0;i44;i++){infileJ1[i];infile.ignore(1);infileI1[i];infile.ignore(1);infileX[i];infile.ignore(1);infileY[i];infile.ignore(1);infileZ[i];infile.ignore(1);}}infile.close();}coutJ1[1]J1[43]endl;//像素坐标转化为直角坐标for(intj=0;j44;j++){x1[j]=(J1[j]-5616/2)*6.410256/1000;z1[j]=(3744/2-I1[j])*6.410256/1000;}coutx1[1]x1[43]endl;//左右相片外方位元素初始值及摄影机主距doubleXs1,Ys1,Zs1,Q1,W1,K1,f1,x01,z01,k11,k21;doubleXs2,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;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;doublerr[N]={0.0},dx[N]={0.0},dz[N]={0.0};doubleXX1,YY1,ZZ1,XX2,YY2,ZZ2;//定义旋转矩阵，系数矩阵,常数项和改正数doubleR1[3][3]={0.0},R2[3][3]={0.0},A1[N][N]={0.0},l1[N][N]={0.0},d[N][N]={0.0};intt=0;coutXs1endl;//组成旋转矩阵R1[0][0]=cos(Q1)*cos(K1)-sin(Q1)*sin(W1)*sin(K1);R1[0][1]=cos(W1)*sin(Q1);R1[0][2]=-cos(Q1)*sin(K1)-sin(Q1)*sin(W1)*cos(K1);R1[1][0]=-sin(Q1)*cos(K1)-cos(Q1)*sin(W1)*sin(K1);R1[1][1]=cos(Q1)*cos(W1);R1[1][2]=sin(Q1)*sin(K1)-cos(Q1)*sin(W1)*cos(K1);R1[2][0]=cos(W1)*sin(K1);R1[2][1]=-sin(W1);R1[2][2]=cos(W1)*cos(K1);R2[0][0]=cos(Q2)*cos(K2)-sin(Q2)*sin(W2)*sin(