近景摄影测量光束法平差报告.docx
《近景摄影测量光束法平差报告.docx》由会员分享，可在线阅读，更多相关《近景摄影测量光束法平差报告.docx（18页珍藏版）》请在冰豆网上搜索。
近景摄影测量光束法平差报告
近景摄影测量
光束法平差报告
2011年6月4日
1作业目的------------------------------------------------------------------------------------3
2外业控制点的观测与解算--------------------------------------------------------3
3近景影像获取----------------------------------------------------------------------------4
4LPS刺点点位------------------------------------------------------------------4
5光束法平差与精度评定------------------------------------------------------------5
6总结---------------------------------------------------------------------------------------------11
1作业目的
以近景摄影测量大实习为基础，对所摄取近景相片解析处理，以外业控制点的解算成果以及内业LPS平差结果为依据，编写光束法平差程序，由22个控制点的像素坐标及5个“已知控制点”的三维坐标求解其余17个控制点的三维坐标，并评定精度。
2作业条件及数据
点号
像素坐标列（J）
像素坐标行（I）
X
Y
Z
左片：
2
650.989
2114.93
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.582
0
0
0
3
660.452
1441.411
0
0
0
4
728.563
816.585
0
0
0
5
1965.895
2557.996
0
0
0
6
1910.105
1210.07
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.7473
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.882
0
0
0
3
685.403
1472.574
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阵为像点坐标误差改正，88行*1列；
A阵为系数阵，88行*73列；
X阵伟未知数阵，73行*1列（未知数包括12个外方位元素、6个内方位元素、51个加密点三维坐标、4个畸变参数）
L阵为常数项阵，88行*1列
由式X=（AtA）-1ATL解求未知数即可
4程序
#include
#include
#include
#include
#include
constintN=90;
//求转置矩阵
templatevoidTranspose（T1*mat1,T2*mat2,inta,intb）
{inti,j;
for（i=0;ifor（j=0;jmat2[j][i]=mat1[i][j];
return;
}
//求矩阵的乘积
templatevoidArray_mul（T1*mat1,T2*mat2,T2*result,inta,intb,intc）
{inti,j,k;
for（i=0;i{for（j=0;j{result[i][j]=0;
for（k=0;kresult[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;i<2*m;i++）
for（j=0;j<2*m;j++）
{if（i==j）C[i][j]=1.0;
elseC[i][j]=0.0;}
for（i=0;ifor（j=0;jB[i][j]=A[i][j];
for（i=0;ifor（j=m;j<2*m;j++）
B[i][j]=C[i][j-m];
cout.precision（5）;
for（k=0;k{for（i=k;ifor（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;kfor（i=0;ifor（j=2*m-1;j>=i;j--）
B[i][j]=B[i][j]-B[i][k]*B[k][j];
for（i=0;ifor（j=0;jA[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;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（intj=0;j<44;j++）
{
x1[j]=（J1[j]-5616/2）*6.410256/1000;
z1[j]=（3744/2-I1[j]）*6.410256/1000;
}
cout<
//左右相片外方位元素初始值及摄影机主距
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;
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（intk=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（intp=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]-x01））;
dz[p]=（z1[p]-z01）*（k11*（z1[p]-z01）*（z1[p]-z01）+k21*（z1[p]-z01）*（z1[p]-z01）*（z1[p]-z01）*（z1[p]-z01））;
}
cout<for（intq=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]-x02））;
dz[q]=（z1[q]-z02）*（k12*（z1[q]-z02）*（z1[q]-z02）+k22*（z1[q]-z02）*（z1[q]-z02）*（z1[q]-z02）*（z1[q]-z02））;
}
cout<//计算系数阵（88*73）和常数项
//左片
for（inti=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*cos（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-f1*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
升级会员 