近景摄影测量光束法平差报告
近景摄影测量光束法平差报告
2011 年 6 月 4 日
1 作业目的------------------------------------------------------------------------------------ 3
2 外业控制点的观测与解算-------------------------------------------------------- 3
3 近景影像获取---------------------------------------------------------------------------- 4
4 LPS刺点点位------------------------------------------------------------------ 4
5 光束法平差与精度评定------------------------------------------------------------ 5
6 总结--------------------------------------------------------------------------------------------- 11
1 作业目的
以近景摄影测量大实习为基础，对所摄取近景相片解析处理，以外业控制点的解算成果以及内业LPS平差结果为依据，编写光束法平差程序，由22个控制点的像素坐标及5个“已知控制点”的三维坐标求解其余17个控制点的三维坐标，并评定精度。
2 作业条件及数据
点号像素坐标列（J）像素坐标行（I）X Y Z
左片：
2 650.989 2114.9
3 497.4532 353.7473 299.8953
8 2792.491 2259.531 508.8008 342.3524 298.6832
10 2791.483 740.514 508.8138 342.3548 307.0717
16 3928.559 2120.49 520.2969 353.7531 300.1146
21 4890.584 2130.45 527.9857 353.5821 300.1037
1 648.624 2765.58
2 0 0 0
3 660.452 1441.411 0 0 0
4 728.563 816.58
5 0 0 0
5 1965.895 2557.99
6 0 0 0
6 1910.105 1210.0
7 0 0 0
7 2767.455 3044.531 0 0 0
9 2774.059 1493.061 0 0 0
12 3319.011 2665.417 0 0 0
13 3312.286 1986.582 0 0 0
14 3298.468 1284.901 0 0 0
15 4055.052 2705.029 0 0 0
17 3808.985 1539.018 0 0 0
18 3715.006 962.032 0 0 0
19 3836.444 706.426 0 0 0
20 4883.39 2691.651 0 0 0
22 4754 1681 0 0 0
23 4825.409 1018.545 0 0 0
右片：
2 670.948 2129.967 497.4532 353.747
3 299.8953
8 2346.443 2264.542 508.8008 342.3524 298.6832
10 2361.448 691.079 508.8138 342.3548 307.0717
16 4088.419 2115.427 520.2969 353.7531 300.1146
20 5203.441 2736.112 527.9857 353.5821 300.1037
1 666.103 2764.88
2 0 0 0
3 685.403 1472.57
4 0 0 0
4 754.414 860.656 0 0 0
5 1652.431 2568.503 0 0 0
6 1600.058 1207.014 0 0 0
7 2312.027 3077.964 0 0 0
9 2334.472 1470.473 0 0 0
12 3083.193 2691.257 0 0 0
13 3083.066 1976.987 0 0 0
14 3072.428 1240.419 0 0 0
15 4230.527 2741.493 0 0 0
17 3956.445 1498.102 0 0 0
18 3852.033 888.02 0 0 0
19 4075.943 613.315 0 0 0
21 5214.492 2119.571 0 0 0
22 5144.463 1507.538 0 0 0
23 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, 相片主距：50mm
3 平差思想
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 L
V阵为像点坐标误差改正，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
#include
#include
#include
#include
const int N=90;
//求转置矩阵
template
{int i,j;
for(i=0;i for(j=0;j mat2[j][i]=mat1[i][j]; return; } //求矩阵的乘积 template for(i=0;i {for(j=0;j {result[i][j]=0; for(k=0;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 for(j=0;j B[i][j]=A[i][j]; for(i=0;i for(j=m;j<2*m;j++) B[i][j]=C[i][j-m]; cout.precision(5); for(k=0;k {for(i=k;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 for(i=0;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 for(j=0;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< //像素坐标转化为直角坐标 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< //左右相片外方位元素初始值及摄影机主距 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< //组成旋转矩阵 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< 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< //计算系数阵(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/X A1[H][1]=(R1[1][1]*(x1[i]-x01)-R1[1][0]*f1)/YY1;//x/Ys A1[H][2]=(R1[2][1]*(x1[i]-x01)-R1[2][0]*f1)/YY1;//x/Zs A1[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/Q A1[H][4]=-f1*sin(K1)-(x1[i]-x01)*((x1[i]-x01)*sin(K1)+(z1[i]-z01)*cos(K1))/f1;//x/W A1[H][5]=z1[i]-z01;//x/K A1[H][6]=(x1[i]-x01)/f1;//x/f A1[H][7]=1;//x/x0 A1[H][8]=0;//x/z0 A1[H][9]=(x1[i]-x01)*rr[i];//x/k1 A1[H][10]=(x1[i]-x01)*rr[i]*rr[i];//x/k2 A1[H][11]=(R2[0][1]*(x1[i]-x02)-R2[0][0]*f2)/YY2;//x/Xs=-x/X A1[H][12]=(R2[1][1]*(x1[i]-x02)-R2[1][0]*f2)/YY2;//x/Ys A1[H][13]=(R2[2][1]*(x1[i]-x02)-R2[2][0]*f2)/YY2;//x/Zs A1[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/Q A1[H][15]=-f2*sin(K2)-(x1[i]-x02)*((x1[i]-x02)*sin(K2)+(z1[i]-z02)*cos(K2))/f2;//x/W A1[H][16]=z1[i]-z02;//x/K A1[H][17]=(x1[i]-x02)/f2;//x/f A1[H][18]=1;//x/x0 A1[H][19]=0;//x/z0 A1[H][20]=(x1[i]-x02)*rr[i];//x/k1 A1[H][21]=(x1[i]-x02)*rr[i]*rr[i];//x/k2 A1[H+1][0]=(R1[0][1]*(z1[i]-z01)-R1[0][2]*f1)/YY1;//z/Xs=-z/X A1[H+1][1]=(R1[1][1]*(z1[i]-z01)-R1[1][2]*f1)/YY1;//z/Ys=-z/Y A1[H+1][2]=(R1[2][1]*(z1[i]-z01)-R1[2][2]*f1)/YY1;//z/Zs=-z/Z A1[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/Q A1[H+1][4]=-f1*cos(K1)-(z1[i]-z01)*((x1[i]-x01)*sin(K1)+(z1[i]-z01)*cos(K1))/f1;//z/W A1[H+1][5]=-(x1[i]-x01);//z/K A1[H+1][6]=(z1[i]-z01)/f1;//z/f A1[H+1][7]=0;//z/x0 A1[H+1][8]=1;//z/z0 A1[H+1][9]=(z1[i]-z01)*rr[i];//x/k1 A1[H+1][10]=(z1[i]-z01)*rr[i]*rr[i];//x/k2 A1[H+1][11]=(R2[0][1]*(z1[i]-z02)-R2[0][2]*f2)/YY2;//z/Xs=-z/X A1[H+1][12]=(R2[1][1]*(z1[i]-z02)-R2[1][2]*f2)/YY2;//z/Ys=-z/Y A1[H+1][13]=(R2[2][1]*(z1[i]-z02)-R2[2][2]*f2)/YY2;//z/Zs=-z/Z A1[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/Q A1[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/K A1[H+1][17]=(z1[i]-z02)/f2;//z/f A1[H+1][18]=0;//z/x0 A1[H+1][19]=1;//z/z0 A1[H+1][20]=(z1[i]-z02)*rr[i];//x/k1 A1[H+1][21]=(z1[i]-z02)*rr[i]*rr[i];//x/k2 H=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/X A1[HH][1]=(R1[1][1]*(x1[ii]-x01)-R1[1][0]*f1)/YY1;//x/Ys A1[HH][2]=(R1[2][1]*(x1[ii]-x01)-R1[2][0]*f1)/YY1;//x/Zs A1[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/Q A1[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/K A1[HH][6]=(x1[ii]-x01)/f1;//x/f A1[HH][7]=1;//x/x0 A1[HH][8]=0;//x/z0 A1[HH][9]=(x1[ii]-x01)*rr[ii];//x/k1 A1[HH][10]=(x1[ii]-x01)*rr[ii]*rr[ii];//x/k2 A1[HH][11]=(R2[0][1]*(x1[ii]-x02)-R2[0][0]*f2)/YY2;//x/Xs=-x/X A1[HH][12]=(R2[1][1]*(x1[ii]-x02)-R2[1][0]*f2)/YY2;//x/Ys A1[HH][13]=(R2[2][1]*(x1[ii]-x02)-R2[2][0]*f2)/YY2;//x/Zs A1[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/Q A1[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/K A1[HH][17]=(x1[ii]-x02)/f2;//x/f A1[HH][18]=1;//x/x0 A1[HH][19]=0;//x/z0 A1[HH][20]=(x1[ii]-x02)*rr[ii];//x/k1 A1[HH][21]=(x1[ii]-x02)*rr[ii]*rr[ii];//x/k2 A1[HH+1][0]=(R1[0][1]*(z1[ii]-z01)-R1[0][2]*f1)/YY1;//z/Xs=-z/X A1[HH+1][1]=(R1[1][1]*(z1[ii]-z01)-R1[1][2]*f1)/YY1;//z/Ys=-z/Y A1[HH+1][2]=(R1[2][1]*(z1[ii]-z01)-R1[2][2]*f1)/YY1;//z/Zs=-z/Z A1[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/Q A1[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/K A1[HH+1][6]=(z1[ii]-z01)/f1;//z/f A1[HH+1][7]=0;//z/x0 A1[HH+1][8]=1;//z/z0 A1[HH+1][9]=(z1[ii]-z01)*rr[ii];//x/k1 A1[HH+1][10]=(z1[ii]-z01)*rr[ii]*rr[ii];//x/k2 A1[HH+1][11]=(R2[0][1]*(z1[ii]-z02)-R2[0][2]*f2)/YY2;//z/Xs=-z/X A1[HH+1][12]=(R2[1][1]*(z1[ii]-z02)-R2[1][2]*f2)/YY2;//z/Ys=-z/Y A1[HH+1][13]=(R2[2][1]*(z1[ii]-z02)-R2[2][2]*f2)/YY2;//z/Zs=-z/Z A1[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/Q A1[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/K A1[HH+1][17]=(z1[ii]-z02)/f2;//z/f A1[HH+1][18]=0;//z/x0 A1[HH+1][19]=1;//z/z0 A1[HH+1][20]=(z1[ii]-z02)*rr[ii];//x/k1 A1[HH+1][21]=(z1[ii]-z02)*rr[ii]*rr[ii];//x/k2 HH=HH+2; } //计算左片外方位元素Xs1,Ys1,Zs1,Q,W,K double 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<<"迭代次数："< cout<<"外方位元素的直线元素为"< cout<<" Xs1 = "<