摄影测量学单像空间后方交会编程实习报告材料Word格式.docx
《摄影测量学单像空间后方交会编程实习报告材料Word格式.docx》由会员分享，可在线阅读，更多相关《摄影测量学单像空间后方交会编程实习报告材料Word格式.docx（10页珍藏版）》请在冰豆网上搜索。
地面点坐标
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"
//矩阵运算头文件，来自网络，已上传
voidmain（）
{
doubleXs,Ys,Zs,q,w,k;
//外方位元素
doublex0,y0,f;
//方位元素
doublex[4],y[4];
//影像坐标
doubleX0[4],Y0[4],Z0[4];
//地面坐标
doublem;
//比例尺
doublea1,a2,a3,b1,b2,b3,c1,c2,c3;
//旋转矩阵
intn=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（inti=0;
i<
4;
i++）
{
Xs+=X0[i];
Ys+=Y0[i];
}
Xs/=4;
Ys/=4;
Zs=f*m;
MatrixA（8,6）;
MatrixL（8,1）;
MatrixAT（6,8）;
MatrixATA（6,6）;
MatrixATA_（6,6）;
MatrixXx（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（inti=0;
i<
4;
i++）
{
doubleX=a1*（X0[i]-Xs）+b1*（Y0[i]-Ys）+c1*（Z0[i]-Zs）;
doubleY=a2*（X0[i]-Xs）+b2*（Y0[i]-Ys）+c2*（Z0[i]-Zs）;
doubleZ=a3*（X0[i]-Xs）+b3*（Y0[i]-Ys）+c3*（Z0[i]-Zs）;
doublexx=x[i]-x0;
doubleyy=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）;
%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）;
//求中误差
doublem0,mi[6];
//单位权中误差的绝对值m0以及第i个未知数的中误差mi
m0=sqrt（（（A*Xx-L）.Transpose（）*（A*Xx-L））.ToDouble（）/（2*4-6））;
MatrixQ（6,6）;
Q=（A.Transpose（）*A）.Inverse（）;
for（inti=0;
6;
mi[i]=（sqrt（Q[i][i]））*m0;
//printf（"
%lf\n"
m0）;
//输出
printf（"
迭代次数:
旋转矩阵R:
\n"
）;
%.5lf%.5lf%.5lf\n"
a1,a2,a3）;
b1,b2,b3）;
c1,c2,c3）;
\n外方位元素解:
\nWs=%.2lfYs=%.2lfZs=%.2lf\nq=%.5lfw=%.5lfk=%.5lf\n\n"
单位权中误差的绝对值:
%lfm\n"
Xs的精度:
mi[0]）;
Ys的精度:
mi[1]）;
Zs的精度:
mi[2]）;
q的精度:
mi[3]）;
w的精度:
mi[4]）;
k的精度:
mi[5]）;
//保存到文件，结果文件默认保存在程序目录
FILE*fp;
fp=fopen（"
result.txt"
"
w"
//fprintf（fp,"
fprintf（fp,"
}
七、运算结果
1.运行结果
2.文件保存结果
八、实习心得
升级会员 