后方交会求外方位元素.docx
《后方交会求外方位元素.docx》由会员分享，可在线阅读，更多相关《后方交会求外方位元素.docx（14页珍藏版）》请在冰豆网上搜索。
后方交会求外方位元素
西南交通大学
摄影测量学课程设计报告
单像空间后方交会
求解外方位元素
目录：
作业任务…………………………………………………………3
计算原理…………………………………………………………3
算法流程…………………………………………………………6
源程序……………………………………………………………7
计算结果…………………………………………………………12
结果分析…………………………………………………………13
心得体会…………………………………………………………13
作业任务：
已知条件
摄影机主距f=153.24mm，x0=0，y0=0,像片比例尺为1:
40000，有四对点的像点坐标与相应的地面坐标如下表。
点号
像点坐标
地面坐标
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
以单像空间后方交会方法，求解该像片的外方位元素。
计算原理：
如图所示，物点A和摄影中心S在地面摄影测量坐标系中的坐标依次是（X,Y,Z）、（XS,YS,ZS）；像点a在像空间坐标系中的坐标是（x,y,-f）。
那么由共线条件方程知，
其中ai，bi，ci是只含三个独立参数ψ，ω，κ的九个方向余弦。
在方程中共有六个未知参数XS,YS,ZS,ψ,ω,κ，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。
由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有
用新的符号表示各偏导数后为
其中（x）、（y）是函数近似值，
是外方位元素近似值的改正数，它们的系数为函数的偏导数。
为了便于推导，令
=
=
=
那么有
对于系数，其严密算法（以
为例）如下：
对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为
为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。
此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数
，可列出每个点的误差方程式
用矩阵表示为
那么由
算法流程：
Ø获取已知数据m,x0,y0,f,Xtp,Ytp,Ztp
Ø量测控制点像点坐标x,y
Ø确定未知数初值Xs0,Ys0,Zs0,0,0,0
Ø组成误差方程式并法化
Ø解求外方位元素改正数
Ø检查迭代是否收敛
具体的流程程序框图如下：
是
否否
是
否
是
源程序：
程序代码：
#include
#include
usingnamespacestd;
constintN=4;
constintn=6;
/*---------------矩阵相乘---------------------*/
voidmult（double*m1,double*m2,double*result,inti_1,intj_12,intj_2）
{
inti,j,k;
for（i=0;ifor（j=0;j{
result[i*j_2+j]=0.0;
for（k=0;kresult[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];
}
return;
}
/*---------------矩阵求逆---------------------*/
voidinverse（doublec[n][n]）
{inti,j,h,k;
doublep;
doubleq[n][12];
for（i=0;ifor（j=0;jq[i][j]=c[i][j];
for（i=0;ifor（j=n;j<12;j++）
{if（i+6==j）
q[i][j]=1;
else
q[i][j]=0;}
for（h=k=0;kfor（i=k+1;i{if（q[i][h]==0）
continue;
p=q[k][h]/q[i][h];
for（j=0;j<12;j++）
{q[i][j]*=p;
q[i][j]-=q[k][j];
}
}
for（h=k=n-1;k>0;k--,h--）//消去对角线以上的数据
for（i=k-1;i>=0;i--）
{if（q[i][h]==0）
continue;
p=q[k][h]/q[i][h];
for（j=0;j<12;j++）
{q[i][j]*=p;
q[i][j]-=q[k][j];}}
for（i=0;i{p=1.0/q[i][i];
for（j=0;j<12;j++）
q[i][j]*=p;}
for（i=0;ifor（j=0;jc[i][j]=q[i][j+6];
}
/*---------------矩阵转置---------------------*/
voidtranspose（double*m1,double*m2,intm,intn）//矩阵转置
{
inti,j;
for（i=0;ifor（j=0;jm2[j*m+i]=m1[i*n+j];
return;
}
voidmain（）
{
doubleXs,Ys,Zs,q,w,k;
doublea[3],b[3],c[3];
doublex0,y0,f;
doublex[N],y[N];
doubleX[N],Y[N],Z[N];
doublex1[N],y1[N];
doublem;
doubleL[2*N];
doubleXX[6];
doubleA[2*N][6];
doubleX0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];
inti,n=0;
doublesum=0,m0;
/*---------------输入点地面坐标---------------------*/
X[0]=36589.41;X[1]=37631.08;X[2]=39100.97;X[3]=40426.54;
Y[0]=25273.32;Y[1]=31324.51;Y[2]=24934.98;Y[3]=30319.81;
Z[0]=2195.17;Z[1]=728.69;Z[2]=2386.50;Z[3]=757.31;
/*---------------输入点像片坐标---------------------*/
x[0]=-86.15;x[1]=-53.40;x[2]=-14.78;x[3]=10.46;
y[0]=-68.99;y[1]=82.21;y[2]=-76.63;y[3]=64.43;
cout</*-----------------设定外方位元素初始值--------------*/
x0=0;y0=0;f=153.24;m=40000;
Xs=0;Ys=0;Zs=f*m/1000;
q=0;w=0;k=0;
XX[3]=1;
/*------------------迭代计算--------------------------*/
while（（XX[3]>6/206265||XX[4]>6/206265||XX[5]>6/206265）&&n<100）
{
/*----------------旋转矩阵R-----------------------*/
a[0]=cos（q）*cos（k）-sin（q）*sin（w）*sin（k）;
a[1]=-cos（q）*sin（k）-sin（q）*sin（w）*cos（k）;
a[2]=-sin（q）*cos（w）;
b[0]=cos（w）*sin（k）;
b[1]=cos（w）*cos（k）;
b[2]=-sin（w）;
c[0]=sin（q）*cos（k）+cos（q）*sin（w）*sin（k）;
c[1]=-sin（q）*sin（k）+cos（q）*sin（w）*cos（k）;
c[2]=cos（q）*cos（w）;
/*-----------------像点坐标计算值------------------*/
for（i=0;i{
X0[i]=a[0]*（X[i]-Xs）+b[0]*（Y[i]-Ys）+c[0]*（Z[i]-Zs）;
Y0[i]=a[1]*（X[i]-Xs）+b[1]*（Y[i]-Ys）+c[1]*（Z[i]-Zs）;
Z0[i]=a[2]*（X[i]-Xs）+b[2]*（Y[i]-Ys）+c[2]*（Z[i]-Zs）;
x1[i]=x0-f*X0[i]/Z0[i];
y1[i]=y0-f*Y0[i]/Z0[i];
}
/*-------------误差方程中各偏导数的值--------------*/
for（i=0;i{
A[2*i][0]=（（a[0]*f+a[2]*（x[i]-x0）））/Z0[i];
A[2*i][1]=（（b[0]*f+b[2]*（x[i]-x0）））/Z0[i];
A[2*i][2]=（（c[0]*f+c[2]*（x[i]-x0）））/Z0[i];
A[2*i][3]=（y[i]-y0）*sin（w）-（（x[i]-x0）*（（x[i]-x0）*cos（k）-y[i]*sin（k））/f+f*cos（k））
*cos（w）;
A[2*i][4]=-f*sin（k）-（x[i]-x0）*（（x[i]-x0）*sin（k）+（y[i]-y0）*cos（k））/f;
A[2*i][5]=y[i]-y0;
L[2*i]=x[i]-x1[i];
A[1+2*i][0]=（（a[1]*f+a[2]*（y[i]-y0）））/Z0[i];
A[1+2*i][1]=（（b[1]*f+b[2]*（y[i]-y0）））/Z0[i];
A[1+2*i][2]=（（c[1]*f+c[2]*（y[i]-y0）））/Z0[i];
A[1+2*i][3]=-（x[i]-x0）*sin（w）-（（y[i]-y0）*（（x[i]-x0）*cos（k）-（y[i]-y0）*sin（k））/f-f*sin（k））
*cos（w）;
A[1+2*i][4]=-f*cos（k）-（y[i]-y0）*（（x[i]-x0）*sin（k）+（y[i]-y0）*cos（k））/f;
A[1+2*i][5]=-x[i]+x0;
L[1+2*i]=y[i]-y1[i];
}
/*-------------------解法方程--------------------*/
transpose（&A[0][0],&At[0][0],2*N,6）;
mult（&At[0][0],&A[0][0],&result1[0][0],6,2*N,6）;
inverse（result1）;
mult（&At[0][0],L,&result2[0][0],6,2*N,1）;
mult（&result1[0][0],&result2[0][0],&XX[0],6,6,1）;
Xs+=XX[0];
Ys+=XX[1];
Zs+=XX[2];
q+=XX[3];
w+=XX[4];
k+=XX[5];
n++;
}
/*----------------旋转矩阵R-----------------------*/
cout<<"迭代次数为:
"<printf（"\n像片外方位元素的解\n"）;
cout<<"航向顷角q:
"<cout<<"旁向倾角w:
"<cout<<"像片旋角k:
"<cout<<"Xs"<cout<<"Ys"<cout<<"Zs"<cout<cout<<"所有坐标已经内置,所以不显示了"<计算结果:
XS=39795.5
YS=27476.5
ZS=7572.69
ψ=-0.00398694
ω=0.00211388
κ=-0.067578
结果分析：
由计算结果可知，最后的计算结果是符合6秒的限差要求的。
由于其中误差未计算，所以精度指标还不明确。
心得体会：
这次的单张航片后方交会计算外方位元素的实习让我对后方交会的原理和解算过程有了更深入的了解，同时对与之相关的一系列知识进行了一个回顾复习。
升级会员 