近景摄影测量光束法平差报告2011 年 6 月 4 日1 作业目的------------------------------------------------------------------------------------ 32 外业控制点的观测与解算-------------------------------------------------------- 33 近景影像获取---------------------------------------------------------------------------- 44 LPS刺点点位------------------------------------------------------------------ 45 光束法平差与精度评定------------------------------------------------------------ 56 总结--------------------------------------------------------------------------------------------- 111 作业目的以近景摄影测量大实习为基础，对所摄取近景相片解析处理，以外业控制点的解算成果以及内业LPS平差结果为依据，编写光束法平差程序，由22个控制点的像素坐标及5个“已知控制点”的三维坐标求解其余17个控制点的三维坐标，并评定精度。

2 作业条件及数据点号像素坐标列（J）像素坐标行（I）X Y Z左片：2 650.989 2114.93 497.4532 353.7473 299.89538 2792.491 2259.531 508.8008 342.3524 298.683210 2791.483 740.514 508.8138 342.3548 307.071716 3928.559 2120.49 520.2969 353.7531 300.114621 4890.584 2130.45 527.9857 353.5821 300.10371 648.624 2765.582 0 0 03 660.452 1441.411 0 0 04 728.563 816.585 0 0 05 1965.895 2557.996 0 0 06 1910.105 1210.07 0 0 07 2767.455 3044.531 0 0 09 2774.059 1493.061 0 0 012 3319.011 2665.417 0 0 013 3312.286 1986.582 0 0 014 3298.468 1284.901 0 0 015 4055.052 2705.029 0 0 017 3808.985 1539.018 0 0 018 3715.006 962.032 0 0 019 3836.444 706.426 0 0 020 4883.39 2691.651 0 0 022 4754 1681 0 0 023 4825.409 1018.545 0 0 0右片：2 670.948 2129.967 497.4532 353.7473 299.89538 2346.443 2264.542 508.8008 342.3524 298.683210 2361.448 691.079 508.8138 342.3548 307.071716 4088.419 2115.427 520.2969 353.7531 300.114620 5203.441 2736.112 527.9857 353.5821 300.10371 666.103 2764.882 0 0 03 685.403 1472.574 0 0 04 754.414 860.656 0 0 05 1652.431 2568.503 0 0 06 1600.058 1207.014 0 0 07 2312.027 3077.964 0 0 09 2334.472 1470.473 0 0 012 3083.193 2691.257 0 0 013 3083.066 1976.987 0 0 014 3072.428 1240.419 0 0 015 4230.527 2741.493 0 0 017 3956.445 1498.102 0 0 018 3852.033 888.02 0 0 019 4075.943 613.315 0 0 021 5214.492 2119.571 0 0 022 5144.463 1507.538 0 0 023 5139.982 903.989 0 0 0外方位元素初始值：(左)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个像点坐标误差方程式组成方程组：V AX LV阵为像点坐标误差改正，88行*1列；A阵为系数阵，88行*73列；X阵伟未知数阵，73行*1列（未知数包括12个外方位元素、6个内方位元素、51个加密点三维坐标、4个畸变参数）L阵为常数项阵，88行*1列由式X=（A t A）-1 A T L解求未知数即可4 程序#include<iomanip.h>#include<stdlib.h>#include<math.h>#include<fstream.h>#include<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++)mat2[j][i]=mat1[i][j];return;}//求矩阵的乘积template<typename T1,typename T2>void Array_mul(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) { int i,j,k;for(i=0;i<a;i++){for(j=0;j<c;j++){result[i][j]=0;for(k=0;k<b;k++)result[i][j]+=mat1[i][k]*mat2[k][j];}}return;}//求逆矩阵inverse(double A[N][N],int m){int i=0,j=0,k=0;double C[20][20],B[20][20];for(i=0;i<2*m;i++)for(j=0;j<2*m;j++){if(i==j) C[i][j]=1.0;else C[i][j]=0.0;}for(i=0;i<m;i++)for(j=0;j<m;j++)B[i][j]=A[i][j];for(i=0;i<m;i++)for(j=m;j<2*m;j++)B[i][j]=C[i][j-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(B[i][k]!=0)B[i][j]=B[i][j]/B[i][k];}for(i=m-1;i>k;i--){if(B[i][k]==0)continue;for(j=0;j<2*m;j++)B[i][j]=B[i][j]-B[k][j];}}for(k=1;k<m;k++)for(i=0;i<k;i++)for(j=2*m-1;j>=i;j--)B[i][j]=B[i][j]-B[i][k]*B[k][j];for(i=0;i<m;i++)for(j=0;j<m;j++)A[j][i]=B[j][m+i];return 1;}void main(){//定义两张相片共个控制点和加密点的像素坐标（J1,I1,J2,I2）和像平面直角坐标（x1,y1,x2,y2），及地面坐标（X,Y,Z）double J1[N]={0.0},I1[N]={0.0},J2[N]={0.0},I2[N]={0.0};double x1[N]={0.0},x2[N]={0.0},z1[N]={0.0},z2[N]={0.0};double X[N]={0.0},Y[N]={0.0},Z[N]={0.0};//导入控制点坐标数据ifstream infile;infile.open("控制点坐标数据.txt");if(infile.is_open()){while(!infile.eof ()){for(int i=0;i<44;i++){infile>>J1[i];infile.ignore(1);infile>>I1[i];infile.ignore(1);infile>>X[i];infile.ignore(1);infile>>Y[i];infile.ignore(1);infile>>Z[i];infile.ignore(1);}}infile.close();}cout<<J1[1]<<" "<<J1[43]<<endl;//像素坐标转化为直角坐标for (int j = 0; j < 44; j++){x1[j] = (J1[j] - 5616/2) * 6.410256 / 1000 ;z1[j] = (3744/2 - I1[j]) * 6.410256 / 1000 ;}cout<<x1[1]<<" "<<x1[43]<<endl;//左右相片外方位元素初始值及摄影机主距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;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 rr[N]={0.0},dx[N]={0.0},dz[N]={0.0};double XX1,YY1,ZZ1,XX2,YY2,ZZ2;//定义旋转矩阵，系数矩阵,常数项和改正数double R1[3][3]={0.0},R2[3][3]={0.0},A1[N][N]={0.0},l1[N][N]={0.0},d[N][N]={0.0};int t=0;cout<<Xs1<<endl;//组成旋转矩阵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(K2);R2[0][1]=cos(W2)*sin(Q2);R2[0][2]=-cos(Q2)*sin(K2)-sin(Q2)*sin(W2)*cos(K2);R2[1][0]=-sin(Q2)*cos(K2)-cos(Q2)*sin(W2)*sin(K2);R2[1][1]=cos(Q2)*cos(W2);R2[1][2]=sin(Q2)*sin(K2)-cos(Q2)*sin(W2)*cos(K2);R2[2][0]=cos(W2)*sin(K2);R2[2][1]=-sin(W2);R2[2][2]=cos(W2)*cos(K2);for(int k=0;k<4;k++){XX1=R1[0][0]*(X[k]-Xs1)+R1[1][0]*(Y[k]-Ys1)+R1[2][0]*(Z[k]-Zs1);YY1=R1[0][1]*(X[k]-Xs1)+R1[1][1]*(Y[k]-Ys1)+R1[2][1]*(Z[k]-Zs1);ZZ1=R1[0][2]*(X[k]-Xs1)+R1[1][2]*(Y[k]-Ys1)+R1[2][2]*(Z[k]-Zs1);XX2=R2[0][0]*(X[k]-Xs1)+R2[1][0]*(Y[k]-Ys1)+R2[2][0]*(Z[k]-Zs1);YY2=R2[0][1]*(X[k]-Xs1)+R2[1][1]*(Y[k]-Ys1)+R2[2][1]*(Z[k]-Zs1);ZZ2=R2[0][2]*(X[k]-Xs1)+R2[1][2]*(Y[k]-Ys1)+R2[2][2]*(Z[k]-Zs1);}for(int p=0;p<22;p++){rr[p]=(x1[p]-x01)*(x1[p]-x01)+(z1[p]-z01)*(z1[p]-z01);dx[p]=(x1[p]-x01)*(k11*(x1[p]-x01)*(x1[p]-x01)+k21*(x1[p]-x01)*(x1[p]-x01)*(x1[p]-x01)*(x1[p]-x0 1));dz[p]=(z1[p]-z01)*(k11*(z1[p]-z01)*(z1[p]-z01)+k21*(z1[p]-z01)*(z1[p]-z01)*(z1[p]-z01)*(z1[p]-z0 1));}cout<<rr[1]<<endl;for(int q=22;q<44;q++){rr[q]=(x1[q]-x02)*(x1[q]-x02)+(z1[q]-z02)*(z1[q]-z02);dx[q]=(x1[q]-x02)*(k12*(x1[q]-x02)*(x1[q]-x02)+k22*(x1[q]-x02)*(x1[q]-x02)*(x1[q]-x02)*(x1[q]-x0 2));dz[q]=(z1[q]-z02)*(k12*(z1[q]-z02)*(z1[q]-z02)+k22*(z1[q]-z02)*(z1[q]-z02)*(z1[q]-z02)*(z1[q]-z0 2));}cout<<rr[22]<<endl;//计算系数阵(88*73)和常数项//左片for(int i=0,H=0;i<22;i++){l1[H][0]=x1[i]-x01+f1*XX1/YY1-dx[i];l1[H+1][0]=z1[i]-z01+f1*ZZ1/YY1-dz[i];A1[H][0]=(R1[0][1]*(x1[i]-x01)-R1[0][0]*f1)/YY1;//x/Xs=-x/XA1[H][1]=(R1[1][1]*(x1[i]-x01)-R1[1][0]*f1)/YY1;//x/YsA1[H][2]=(R1[2][1]*(x1[i]-x01)-R1[2][0]*f1)/YY1;//x/ZsA1[H][3]=(z1[i]-z01)*sin(W1)-((x1[i]-x01)*((x1[i]-x01)*cos(K1)-(z1[i]-z01)*sin(K1))/f1+f1*c os(K1))*cos(W1);//x/QA1[H][4]=-f1*sin(K1)-(x1[i]-x01)*((x1[i]-x01)*sin(K1)+(z1[i]-z01)*cos(K1))/f1;//x/WA1[H][5]=z1[i]-z01;//x/KA1[H][6]=(x1[i]-x01)/f1;//x/fA1[H][7]=1;//x/x0A1[H][8]=0;//x/z0A1[H][9]=(x1[i]-x01)*rr[i];//x/k1A1[H][10]=(x1[i]-x01)*rr[i]*rr[i];//x/k2A1[H][11]=(R2[0][1]*(x1[i]-x02)-R2[0][0]*f2)/YY2;//x/Xs=-x/XA1[H][12]=(R2[1][1]*(x1[i]-x02)-R2[1][0]*f2)/YY2;//x/YsA1[H][13]=(R2[2][1]*(x1[i]-x02)-R2[2][0]*f2)/YY2;//x/ZsA1[H][14]=(z1[i]-z02)*sin(W2)-((x1[i]-x02)*((x1[i]-x02)*cos(K2)-(z1[i]-z02)*sin(K2))/f2+f2* cos(K2))*cos(W2);//x/QA1[H][15]=-f2*sin(K2)-(x1[i]-x02)*((x1[i]-x02)*sin(K2)+(z1[i]-z02)*cos(K2))/f2;//x/WA1[H][16]=z1[i]-z02;//x/KA1[H][17]=(x1[i]-x02)/f2;//x/fA1[H][18]=1;//x/x0A1[H][19]=0;//x/z0A1[H][20]=(x1[i]-x02)*rr[i];//x/k1A1[H][21]=(x1[i]-x02)*rr[i]*rr[i];//x/k2A1[H+1][0]=(R1[0][1]*(z1[i]-z01)-R1[0][2]*f1)/YY1;//z/Xs=-z/XA1[H+1][1]=(R1[1][1]*(z1[i]-z01)-R1[1][2]*f1)/YY1;//z/Ys=-z/YA1[H+1][2]=(R1[2][1]*(z1[i]-z01)-R1[2][2]*f1)/YY1;//z/Zs=-z/ZA1[H+1][3]=-(x1[i]-x01)*sin(W1)-((z1[i]-z01)*((x1[i]-x01)*cos(K1)-(z1[i]-z01)*sin(K1))/f1-f 1*sin(K1))*cos(W1);//z/QA1[H+1][4]=-f1*cos(K1)-(z1[i]-z01)*((x1[i]-x01)*sin(K1)+(z1[i]-z01)*cos(K1))/f1;//z/WA1[H+1][5]=-(x1[i]-x01);//z/KA1[H+1][6]=(z1[i]-z01)/f1;//z/fA1[H+1][7]=0;//z/x0A1[H+1][8]=1;//z/z0A1[H+1][9]=(z1[i]-z01)*rr[i];//x/k1A1[H+1][10]=(z1[i]-z01)*rr[i]*rr[i];//x/k2A1[H+1][11]=(R2[0][1]*(z1[i]-z02)-R2[0][2]*f2)/YY2;//z/Xs=-z/XA1[H+1][12]=(R2[1][1]*(z1[i]-z02)-R2[1][2]*f2)/YY2;//z/Ys=-z/YA1[H+1][13]=(R2[2][1]*(z1[i]-z02)-R2[2][2]*f2)/YY2;//z/Zs=-z/ZA1[H+1][14]=-(x1[i]-x02)*sin(W2)-((z1[i]-z02)*((x1[i]-x02)*cos(K2)-(z1[i]-z02)*sin(K2))/f2-f2*sin(K2))*cos(W2);//z/QA1[H+1][15]=-f2*cos(K2)-(z1[i]-z02)*((x1[i]-x02)*sin(K2)+(z1[i]-z02)*cos(K2))/f2;//z/W A1[H+1][16]=-(x1[i]-x02);//z/KA1[H+1][17]=(z1[i]-z02)/f2;//z/fA1[H+1][18]=0;//z/x0A1[H+1][19]=1;//z/z0A1[H+1][20]=(z1[i]-z02)*rr[i];//x/k1A1[H+1][21]=(z1[i]-z02)*rr[i]*rr[i];//x/k2H=H+2;}//右片for(int ii=22,HH=44;ii<44;ii++){l1[HH][0]=x1[ii]-x01+f1*XX1/YY1-dx[ii];l1[HH+1][0]=z1[ii]-z01+f1*ZZ1/YY1-dz[ii];A1[HH][0]=(R1[0][1]*(x1[ii]-x01)-R1[0][0]*f1)/YY1;//x/Xs=-x/XA1[HH][1]=(R1[1][1]*(x1[ii]-x01)-R1[1][0]*f1)/YY1;//x/YsA1[HH][2]=(R1[2][1]*(x1[ii]-x01)-R1[2][0]*f1)/YY1;//x/ZsA1[HH][3]=(z1[ii]-z01)*sin(W1)-((x1[ii]-x01)*((x1[ii]-x01)*cos(K1)-(z1[ii]-z01)*sin(K1))/f1 +f1*cos(K1))*cos(W1);//x/QA1[HH][4]=-f1*sin(K1)-(x1[ii]-x01)*((x1[ii]-x01)*sin(K1)+(z1[ii]-z01)*cos(K1))/f1;//x/W A1[HH][5]=z1[ii]-z01;//x/KA1[HH][6]=(x1[ii]-x01)/f1;//x/fA1[HH][7]=1;//x/x0A1[HH][8]=0;//x/z0A1[HH][9]=(x1[ii]-x01)*rr[ii];//x/k1A1[HH][10]=(x1[ii]-x01)*rr[ii]*rr[ii];//x/k2A1[HH][11]=(R2[0][1]*(x1[ii]-x02)-R2[0][0]*f2)/YY2;//x/Xs=-x/XA1[HH][12]=(R2[1][1]*(x1[ii]-x02)-R2[1][0]*f2)/YY2;//x/YsA1[HH][13]=(R2[2][1]*(x1[ii]-x02)-R2[2][0]*f2)/YY2;//x/ZsA1[HH][14]=(z1[ii]-z02)*sin(W2)-((x1[ii]-x02)*((x1[ii]-x02)*cos(K2)-(z1[ii]-z02)*sin(K2))/f 2+f2*cos(K2))*cos(W2);//x/QA1[HH][15]=-f2*sin(K2)-(x1[ii]-x02)*((x1[ii]-x02)*sin(K2)+(z1[ii]-z02)*cos(K2))/f2;//x/W A1[HH][16]=z1[ii]-z02;//x/KA1[HH][17]=(x1[ii]-x02)/f2;//x/fA1[HH][18]=1;//x/x0A1[HH][19]=0;//x/z0A1[HH][20]=(x1[ii]-x02)*rr[ii];//x/k1A1[HH][21]=(x1[ii]-x02)*rr[ii]*rr[ii];//x/k2A1[HH+1][0]=(R1[0][1]*(z1[ii]-z01)-R1[0][2]*f1)/YY1;//z/Xs=-z/XA1[HH+1][1]=(R1[1][1]*(z1[ii]-z01)-R1[1][2]*f1)/YY1;//z/Ys=-z/YA1[HH+1][2]=(R1[2][1]*(z1[ii]-z01)-R1[2][2]*f1)/YY1;//z/Zs=-z/ZA1[HH+1][3]=-(x1[ii]-x01)*sin(W1)-((z1[ii]-z01)*((x1[ii]-x01)*cos(K1)-(z1[ii]-z01)*sin(K1)) /f1-f1*sin(K1))*cos(W1);//z/QA1[HH+1][4]=-f1*cos(K1)-(z1[ii]-z01)*((x1[ii]-x01)*sin(K1)+(z1[ii]-z01)*cos(K1))/f1;//z/W A1[HH+1][5]=-(x1[ii]-x01);//z/KA1[HH+1][6]=(z1[ii]-z01)/f1;//z/fA1[HH+1][7]=0;//z/x0A1[HH+1][8]=1;//z/z0A1[HH+1][9]=(z1[ii]-z01)*rr[ii];//x/k1A1[HH+1][10]=(z1[ii]-z01)*rr[ii]*rr[ii];//x/k2A1[HH+1][11]=(R2[0][1]*(z1[ii]-z02)-R2[0][2]*f2)/YY2;//z/Xs=-z/XA1[HH+1][12]=(R2[1][1]*(z1[ii]-z02)-R2[1][2]*f2)/YY2;//z/Ys=-z/YA1[HH+1][13]=(R2[2][1]*(z1[ii]-z02)-R2[2][2]*f2)/YY2;//z/Zs=-z/ZA1[HH+1][14]=-(x1[ii]-x02)*sin(W2)-((z1[ii]-z02)*((x1[ii]-x02)*cos(K2)-(z1[ii]-z02)*sin(K2) )/f2-f2*sin(K2))*cos(W2);//z/QA1[HH+1][15]=-f2*cos(K2)-(z1[ii]-z02)*((x1[ii]-x02)*sin(K2)+(z1[ii]-z02)*cos(K2))/f2;//z/W A1[HH+1][16]=-(x1[ii]-x02);//z/KA1[HH+1][17]=(z1[ii]-z02)/f2;//z/fA1[HH+1][18]=0;//z/x0A1[HH+1][19]=1;//z/z0A1[HH+1][20]=(z1[ii]-z02)*rr[ii];//x/k1A1[HH+1][21]=(z1[ii]-z02)*rr[ii]*rr[ii];//x/k2HH=HH+2;}//计算左片外方位元素Xs1,Ys1,Zs1,Q,W,Kdouble A1T[N][N]={0},A1TA1[N][N]={0},A1TL1[N][N]={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+d[0][0]; Ys1=Ys1+d[1][0]; Zs1=Zs1+d[2][0];Q1=Q1+d[3][0]; W1=W1+d[4][0]; K1=K1+d[5][0];f1=f1+d[6][0]; x01=x01+d[7][0]; z01=z01+d[8][0];k11=k11+d[9][0]; k21=k21+d[10][0];Xs2=Xs2+d[11][0]; Ys2=Ys1+d[12][0]; Zs2=Zs2+d[13][0];Q2=Q2+d[14][0]; W2=W2+d[15][0]; K2=K2+d[16][0];f2=f2+d[17][0]; x02=x02+d[18][0]; z02=z02+d[19][0];k12=k11+d[9][0]; k22=k22+d[10][0];cout<<"迭代次数："<<t<<endl;cout<<"外方位元素的直线元素为"<<endl;cout<<" Xs1 = "<<Xs1<<" , Xs2 = "<<Xs2<<endl;cout<<" Ys1 = "<<Ys1<<" , Ys2 = "<<Ys2<<endl;cout<<" Zs1 = "<<Zs1<<" , Zs2 = "<<Zs2<<endl;cout<<" Q1 = "<<Q1<<" , Q2 = "<<Q2<<endl;cout<<" W1 = "<<Q1<<" , W2 = "<<W2<<endl;cout<<" K1 = "<<K1<<" , K2 = "<<K2<<endl;cout<<endl;cout<<" f1 = "<<f1<<" , f2 = "<<f2<<endl;cout<<" x01 = "<<f1<<" , x02 = "<<f2<<endl;cout<<" z01 = "<<f1<<" , z02 = "<<f2<<endl;cout<<endl;cout<<" k11 = "<<k11<<" , k12 = "<<k12<<endl;cout<<" k21 = "<<k21<<" , k22 = "<<k22<<endl;}5 结果在这里要请求老师原谅，很抱歉在编制程序过程中，出现了陷入死循环的问题，一直修改调试也未能解决问题，导致无法计算出结果。