摄影测量学单像空间后方交会编程实习报告.doc

摄影测量学
单像空间后方交会编程
实习报告
班 级: 130x
姓 名: xx
学 号: 2013302590xxx
指导老师: 李 欣
一、 实习目得
通过对提供得数据进行计算,输出像片得外方位元素并评定精度。深入理解单像空间后方交会得思想,体会在有多余观测情况下,用最小二乘平差方法编程实现解求影像外方位元素得过程。通过尝试编程实现加强编程处理问题得能力与对实习内容得理解,通过对实验结果得分析,增强综合运用所学知识解决实际问题得能力。了解摄影测量平差得基本过程,掌握空间后方交会得定义与实现算法。
二、 实习内容
根据学习得单像空间后方交会得知识,用程序设计语言(C++或C语言)编写一个完整得单像空间后方交会程序,通过对提供得数据进行计算,输出像片得外方位元素并评定精度。
三、 实习数据
已知航摄仪得内方位元素:fk=153、24mm,x0=y0=0,摄影比例尺为1:15000;
4个地面控制点得地面坐标及其对应像点得像片坐标:

像片坐标
地面点坐标
x(mm)
y(mm)
X(m)
Y(m)
Z(m)
1
-86、15
-68、99
36589、41
25273、32
2195、17
2
-53、40
82、21
37631、08
31324、51
728、69
3
-14、78
-76、63
39100、97
24934、98
2386、50
4
10、46
64、43
40426、54
30319、81
757、31
四、 实习原理
如果我们知道每幅影像得6个外方位元素,就能确定被摄物体与航摄影像得关系。因此,如何获取影像得外方位元素,一直就是摄影测量工作者所探讨得问题。可采取得方法有:利用雷达、全球定位系统(GPS)、惯性导航系统(INS)以及星相摄影机来获取影像得外方位元素;也可以利用影像覆盖范围内一定数量得控制点得空间坐标与摄影坐标,根据共线条件方程,反求该影像得外方位元素,这种方法称为单幅影像得空间后方交会。 
单像空间后方交会得基本思想就是:以单幅影像为基础,从该影像所覆盖地面范围内若干控制点得已知地面坐标与相应点得像坐标量测值出发,根据共线条件方程,解求该影像在航空摄影时刻得外方位元素Xs,Ys,Zs,ϕ,ω,κ。
五、 实习流程
1. 获取已知数据。从摄影资料中查取影像比例尺1/m,平均摄影距离(航空摄影得航高、内方位元素x0,y0,;获取控制点得空间坐标Xt,Yt,Zt。 
2. 量测控制点得像点坐标并进行必要得影像坐标系统误差改正,得到像点坐标。
3. 确定未知数得初始值。单像空间后方交会必须给出待定参数得初始值,在竖直航空摄影且地面控制点大体对称分布得情况下,可按如下方法确定初始值:

式中:m为摄影比例尺分母,n为控制点个数; 
4. 计算旋转矩阵R。利用角元素得近似值按公式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3,组成R阵。 
5. 逐点计算像点坐标得近似值。利用未知数得近似值按共线条件方程计算控制点像点坐标得近似值(x),(y)。 
6. 逐点计算误差方程式得系数与常数项,组成误差方程。
7. 计算法方程得系数矩阵ATA与常数项ATL,组成法方程。 
8. 解求外方位元素。根据法方程,解求外方位元素得改正数,并与相应得近似值求与,得到外方位元素新得近似值。
9. 检查计算就是否收敛。将所求得得外方位元素得改正数与规定得限差比较,通常对ϕ,ω,κ得改正数给予限差,当改正数小于限差时,迭代结束。否则用新得近似值重复(4)——(8)步骤计算,直到满足要求为止。 
10. 空间后方交会得精度估计: 
按上述方法所求得得影像外方位元素得精度可以通过法方程式中未知数得系数矩阵得逆阵(ATA)-1来解求,此时视像点坐标为等精度不相关观测值。因为(ATA)-1中第i个主对角线上得元素Qii就就是法方程式中第i个未知数得权倒数,若单位权中误差为m0,则第i个未知数得中误差为:

当参加空间后方交会得控制点有n个时,则单位权中误差可按下式计算:

流程图如下:

六、 程序代码
#include <stdio、h>
#include "Matrix、h" //矩阵运算头文件,来自网络,已上传
void main()
{
double Xs,Ys,Zs,q,w,k; //外方位元素
double x0,y0,f; //内方位元素
double x[4],y[4]; //影像坐标
double X0[4],Y0[4],Z0[4]; //地面坐标
double m; //比例尺
double a1,a2,a3,b1,b2,b3,c1,c2,c3; //旋转矩阵
int n=0; //迭代次数
//地面坐标
X0[0]=36589、41;
X0[1]=37631、08;
X0[2]=39100、97;
X0[3]=40426、54;
Y0[0]=25273、32;
Y0[1]=31324、51;
Y0[2]=24934、98;
Y0[3]=30319、81;
Z0[0]=2195、17;
Z0[1]=728、69;
Z0[2]=2386、50;
Z0[3]=757、31;
//影像坐标
x[0]=-0、08615;
x[1]=-0、05340;
x[2]=-0、01478;
x[3]=0、01046;
y[0]=-0、06899;
y[1]=0、08221;
y[2]=-0、07663;
y[3]=0、06443;
//确定内外方位元素初始值
x0=0;y0=0;
f=153、24/1000;
m=15000;
Xs=0;Ys=0;
q=0;w=0;k=0;
for (int i=0;i<4;i++)
{
Xs += X0[i];
Ys += Y0[i];
}
Xs /= 4;
Ys /= 4;
Zs=f*m;
Matrix A(8, 6);
Matrix L(8, 1);
Matrix AT(6, 8);
Matrix ATA(6, 6);
Matrix ATA_(6, 6);
Matrix Xx(6, 1);
//迭代计算
do
{
//旋转矩阵R
a1 = cos(q)*cos(k) - sin(q)*sin(w)*sin(k);
a2 = -cos(q)*sin(k) - sin(q)*sin(w)*cos(k);
a3 = -sin(q)*cos(w);
b1 = cos(w)*sin(k);
b2 = cos(w)*cos(k);
b3 = -sin(w);
c1 = sin(q)*cos(k) + cos(q)*sin(w)*sin(k);
c2 = -sin(q)*sin(k) + cos(q)*sin(w)*cos(k);
c3 = cos(q)*cos(w);
//计算像点坐标
for (int i = 0; i < 4; i++)
{
double X = a1 * ( X0[i] - Xs) + b1 * (Y0[i] - Ys) + c1 * (Z0[i] - Zs);
double Y = a2 * ( X0[i] - Xs) + b2 * (Y0[i] - Ys) + c2 * (Z0[i] - Zs);
double Z = a3 * ( X0[i] - Xs) + b3 * (Y0[i] - Ys) + c3 * (Z0[i] - Zs);
double xx = x[i]-x0;
double yy = y[i]-y0;
A[2 * i][0] = (a1*f + a3*xx) / Z;
A[2 * i][1] = (b1*f + b3*xx) / Z;
A[2 * i][2] = (c1*f + c3*xx) / Z;
A[2 * i][3] = yy*sin(w) - (xx*(xx*cos(k) - yy*sin(k)) / f + f*cos(k))*cos(w);
A[2 * i][4] = -f*sin(k) - xx*(xx*sin(k) + yy*cos(k)) / f;
A[2 * i][5] = yy;
A[2 * i + 1][0] = (a2*f + a3*yy) / Z;
A[2 * i + 1][1] = (b2*f + b3*yy) / Z;
A[2 * i + 1][2] = (c2*f + c3*yy) / Z;
A[2 * i + 1][3] = -xx*sin(w) - (yy*(xx*cos(k) - yy*sin(k)) / f - f*sin(k))*cos(w);
A[2 * i + 1][4] = -f*cos(k) - yy*(xx*sin(k) + yy*cos(k)) / f;
A[2 * i + 1][5] = -xx;
L[2 * i][0] = x[i] - (-f*X / Z);
L[2 * i + 1][0] =y[i] - (-f*Y / Z);
}
//法方程得解X=(ATA)_ATL
AT = A、Transpose(); //转置
ATA = AT * A;
ATA_ = ATA、Inverse(); //求逆
Xx = (ATA_ * AT) * L; //法方程解
Xs += Xx[0][0];
Ys += Xx[1][0];
Zs += Xx[2][0];
q += Xx[3][0];
w += Xx[4][0];
k += Xx[5][0];
n++;
//printf("%d\n",n);
//printf("%lf %lf %lf %、5lf %、5lf %、5lf\n", Xs, Ys, Zs, q, w, k);
}while ((abs(Xx[3][0]) > 0、000001 || abs(Xx[4][0]) > 0、000001 || abs(Xx[5][0]) > 0、000001) && n<100);
//求中误差
double m0,mi[6]; //单位权中误差得绝对值m0以及第i个未知数得中误差mi
m0 = sqrt(((A*Xx - L)、Transpose()*(A*Xx - L))、ToDouble() / (2 * 4 - 6));
Matrix Q(6,6);
Q = (A、Transpose()*A)、Inverse();
for (int i = 0; i < 6; i++)
{
mi[i] = (sqrt(Q[i][i]))*m0;
}
//printf("%lf\n",m0);
//输出
printf("迭代次数:%d\n",n);
printf("旋转矩阵R:\n");
printf("%、5lf %、5lf %、5lf\n",a1,a2,a3);
printf("%、5lf %、5lf %、5lf\n",b1,b2,b3);
printf("%、5lf %、5lf %、5lf\n",c1,c2,c3);
printf("\n外方位元素解:\nWs=%、2lf Ys=%、2lf Zs=%、2lf\nq=%、5lf w=%、5lf k=%、5lf\n\n", Xs, Ys, Zs, q, w, k);
printf("单位权中误差得绝对值:%lfm\n",m0);
printf("Xs得精度:%lfm\n",mi[0]);
printf("Ys得精度:%lfm\n",mi[1]);
printf("Zs得精度:%lfm\n",mi[2]);
printf("q得精度:%lf\n",mi[3]);
printf("w得精度:%lf\n",mi[4]);
printf("k得精度:%lf\n",mi[5]);
//保存到文件,结果文件默认保存在程序目录
FILE* fp;
fp=fopen("result、txt","w");
// fprintf(fp,"迭代次数:%d\n",n);
fprintf(fp,"旋转矩阵R:\n");
fprintf(fp,"%、5lf %、5lf %、5lf\n",a1,a2,a3);
fprintf(fp,"%、5lf %、5lf %、5lf\n",b1,b2,b3);
fprintf(fp,"%、5lf %、5lf %、5lf\n",c1,c2,c3);
fprintf(fp,"\n外方位元素解:\nWs=%、2lf Ys=%、2lf Zs=%、2lf\nq=%、5lf w=%、5lf k=%、5lf\n\n", Xs, Ys, Zs, q, w, k);
fprintf(fp,"单位权中误差得绝对值:%lfm\n",m0);
fprintf(fp,"Xs得精度:%lfm\n",mi[0]);
fprintf(fp,"Ys得精度:%lfm\n",mi[1]);
fprintf(fp,"Zs得精度:%lfm\n",mi[2]);
fprintf(fp,"q得精度:%lf\n",mi[3]);
fprintf(fp,"w得精度:%lf\n",mi[4]);
fprintf(fp,"k得精度:%lf\n",mi[5]);
}
七、 运算结果
1. 运行结果

2. 文件保存结果

八、 实习心得